EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 699 (July 2025) A&A, 699 (2025) A71 Full HTML
Open Access
Issue
A&A
Volume 699, July 2025
Article Number A71
Number of page(s) 17
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/202453472
Published online 01 July 2025

© The Authors 2025

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.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

The mass lost from very massive stars (VMSs) with masses over 100 M⊙ at low metallicity is a key factor that shapes the evolution of these objects and their chemical yields. Unlike stars at solar metallicity, where line-driven winds dominate the mass loss, stars in low-Z environments experience significantly weaker winds (Vink et al. 2001), and this leads to less efficient angular momentum loss and a potentially more prolonged phase of rapid rotation, which can enhance internal mixing and alter nucleosynthesis outcomes. Studies have suggested that VMSs at low metallicity could retain much of their mass until late evolutionary stages, resulting in distinct chemical feedback patterns compared to their higher-Z counterparts (Yusof et al. 2013).

The chemical evolution of globular clusters (GCs) has long been a subject of intense study, particularly due to the well-established abundance anomalies such as the Na-O anti-correlation (Carretta et al. 2009; Charbonnel 2016; Bastian & Lardo 2018). These anomalies suggest that GCs experienced self-enrichment processes early in their formation, likely driven by feedback from a first population of stars. Understanding the sources of this chemical enrichment requires detailed models of stellar nucleosynthesis and mass loss, particularly for stars at low metallicity (Z). The focus on low-Z environments is motivated by their relevance to both the early Universe and nearby dwarf galaxies, which are considered analogues of the conditions in the first galaxies.

The initial mass function (IMF) in these low-metallicity environments is thought to be more top heavy, with a higher proportion of very massive stars than is typically seen in the present-day Universe. This expectation follows from theoretical studies of star formation under primordial conditions (Bromm & Larson 2004) and diminished wind feedback Vink (2018). The first empirical evidence of a more top-heavy IMF, found in the Large Magellanic Cloud (LMC, ∼50% solar Z) was presented in Schneider et al. (2018). The feedback from these VMSs is considered to play a crucial role in driving early galaxy evolution, providing both chemical enrichment and energy to their surroundings (Krause et al. 2013; Vink et al. 2015; Crowther et al. 2016). In dwarf galaxies, the presence of abundance anomalies similar to those found in GCs points to a similar population of massive stars contributing to the chemical evolution (Gratton et al. 2004).

Previous works have focused on the nucleosynthetic yields of massive stars, particularly at low metallicity, and their potential to explain the observed abundance patterns in both GCs and the broader Universe. These studies have emphasised the role of rotational mixing at low Z, which may be able to bring products of nuclear burning to the surface, thereby contributing to the observed Na-O anti-correlation (Decressin et al. 2007). It should be noted, however, that such a model relies not only on high stellar rotation rates but also on very efficient rotational mixing, which has been questioned (Vink & Harries 2017). Denissenkov & Hartwick (2014) proposed supermassive stars (M* > 104 M⊙) as the source for observed anti-correlations in GCs due to losing significant amounts of mass within the first 10% of core H burning. However, Denissenkov & Hartwick (2014) also suggested that super-Eddington winds would be required to provide such high mass loss on this timescale, which crucially occurs before further nucleosynthesis.

Instead, the early Universe was likely populated by a significant number of VMSs with masses in the range 100–1000 M⊙, which may have contributed to the reionisation and chemical enrichment of the intergalactic medium (Schaerer et al. 2025). These stars are thought to have provided the first significant chemical feedback in the Universe, enriching primordial gas and setting the stage for subsequent generations of star formation (Baraffe et al. 2001). In this context, studying the yields of very massive stars at low metallicity not only helps bring greater understanding of the formation of GCs but also provides insights into the broader processes that shaped the early Universe.

This paper aims to investigate the wind yields from very massive stars at low metallicity, focusing on their contributions to chemical feedback in GCs and dwarf galaxies. We explore how different mass-loss rates, rotation speeds, and metallicities influence the chemical enrichment patterns, particularly in the context of Na-O anti-correlations and other abundance anomalies. By comparing our theoretical yields with recent observations, we aim to provide new insights into the role of VMSs in the early chemical evolution of galaxies.

2. Stellar models

We calculate grids of models using the one-dimensional stellar evolution code MESA (r10398; Paxton et al. 2011, 2013, 2015, 2018, 2019) for a grid of initial masses of 30, 50, 80, 100, 200, 300, 400, and 500 M⊙. We implement a nuclear reaction network which includes the relevant isotopes for massive star evolution until the end of core O burning such that the nucleosynthesis is comparable to Higgins et al. (2023). This nuclear network comprises the following 92 isotopes: n, 1, 2H, 3, 4He, 6, 7Li, 7, 9, 10Be, 8, 10, 11B, 12, 13C, 13 − 16N, 14 − 19O, 17 − 20F, 18 − 23Ne, 21 − 24Na, 23 − 27Mg, 25 − 28Al, 27 − 33Si, 30 − 34P, 31 − 37S, 35 − 38Cl, 35 − 41Ar, 39 − 44K, and 39 − 44, 46, 48Ca. Our stellar models are computed for a range of metallicities, ranging from the metallicity of the Small Magellanic Cloud (ZSMC) to 1% Z⊙ where Z⊙ = 0.02, X = 0.720, and Y = 0.26, where the relative composition is adopted from Asplund et al. (2009). At each metallicity we calculate Y = Yprim + (ΔY/ΔZ)Z where ΔY/ΔZ = 2 and Yprim = 0.24 (Anders & Grevesse 1989; Pols et al. 1995), see Table 1. We compare with the Z⊙ models of Higgins et al. (2023) in Figs. 2 and 3, and we note that in that case, models were calculated with Z⊙ = 0.014. We used the OPAL opacity tables from Rogers & Nayfonov (2002) and adopted nuclear reaction rates from the JINA Reaclib Database v2.2 (Cyburt et al. 2010). The mixing-length-theory (MLT) of convection describes the treatment of convection in our models, where we apply an efficiency of αmlt = 1.67 (Arnett et al. 2019). The Schwarzschild criterion defines the convective boundaries in our models and as such we do not implement semiconvective mixing. For convective boundary mixing (CBM), we include the exponential decaying diffusive model of Freytag et al. (1996) (see also Herwig 2000) with fov = 0.03 (corresponding to αov ≃ 0.3) for the top of convective cores and shells, and with fov = 0.006 for the bottom of convective shells. Bowman (2020) find a range of αov from asteroseismology results with αov up to 0.4, so our value for the top of convective zones falls in the range of asteroseismology-inferred values. This value also falls in between the majority of published large grids of stellar models such as αov = 0.1 in Ekström et al. (2012), αov = 0.335 in Brott et al. (2011), and recent studies on CBM (Higgins & Vink 2019; Scott et al. 2021) supporting values for αov up to at least 0.5 for stars above 20 M⊙. For the bottom boundary, a CBM value of 1/5 the value of the top boundary is based on 3D hydrodynamic simulations (Cristini et al. 2017, 2019; Rizzuti et al. 2022) finding that CBM is slower at the bottom boundaries due to being stiffer and therefore harder to penetrate.

Table 1.

Initial composition of our model grid for a range of initial Z in mass fractions.

We apply convection in superadiabatic layers via the MLT++ prescription which aids convergence of such models. The temporal resolution of our models has been set with varcontrol target = 0.0001, and a corresponding spatial resolution of mesh delta = 0.5. All calculations begin with a pre-main sequence (MS) and then evolve from the zero age main sequence (ZAMS) until core H-depletion ( 1 H c < 0.0001 $ ^{1}\rm{H}_\mathrm{{c}} < 0.0001 $) or core He-depletion ( 4 He c < 0.0001 $ ^{4}\rm{He}_\mathrm{{c}} < 0.0001 $). At Z < 10% Z⊙ stellar models evolve at higher temperatures and luminosities, and as such are compact and evolve close to their Eddington limit. Therefore, we calculate the core H-burning MS only for the grids of models with Z < 10% Z⊙ (i.e. Z = 0.00067 and 0.0002), while the higher Z models evolve until core He exhaustion (10% and 20% Z⊙ grids).

The effect of rotation on nucleosynthesis and wind yields are tested for models with Z = 0.004 and 0.002. We focus on non-rotating models at the lowest Z (1/30th Z⊙ and 1/100th Z⊙) since these stars can spin up to reach criticality. However, we provide a subset of models initially rotating at 40% of critical, calculated at 1/30th Z⊙ to compare with observed Na and O abundance patterns (see Sect. 6.1). To explore the consequences of rotation on wind yields at Z = 0.004 and 0.002, we calculate three sets of models for each of these Z: non-rotating, Ωini/Ωcrit = 40%, and 70%. We adopt rotation with instabilities employed for angular momentum transfer and chemical mixing as described by Heger et al. (2000), and do not include the effects of a magnetic field.

thumbnail Fig. 1.

Mass-loss rates of 100 M⊙ models for each Z as a function of core He abundance during the MS evolution.

The wind rates applied throughout our grids of models are invoked for each evolutionary phase and based on the wind driving mechanism at each stage. Therefore, for hot (log10 (Teff/K) ≥ 4.0), optically thin OB stars, evolving comfortably below the transition point (M ≲ 60 M⊙), we adopt the Vink et al. (2001) rates.

We applied the mass-loss framework of Sabhahit et al. (2023) where stars switch from optically thin to optically thick rates above the transition point. This is calculated iteratively for L and M, where the condition for switching to an optically thick wind is established as a function of the wind efficiency parameter η. Above the switch, or transition point, we then apply wind rates based on the relation by Vink et al. (2011).

From 4 kK < (Teff/K) < 100 kK, we employed the same mass-loss framework on the post-MS as during the MS since the core evolution should not dictate the wind driving mechanism. However, below Teff ≈ 4 kK we switch to the cool supergiant prescription with rates adopted from de Jager et al. (1988). Finally, for classical Wolf-Rayet stars, where Teff > 100 kK and Xc < 0.01, we apply the wind relation from the hydrodynamically consistent models of Sander & Vink (2020).

For rotating models, we assume that near the radiative Eddington limit Γ > 0.639, the ΩΓ-limit can be reached, and increased mass can be lost, accounting for the von Zeipel theorem. The implementation of this mass-loss relation follows the approximation as derived by Maeder & Meynet (2000):

M ˙ ( Ω ) M ˙ ( Ω = 0 ) ≈ ( 1 − Γ ) 1 α − 1 ( 1 − 4 9 ( v v crit , 1 ) 2 − Γ ) 1 α − 1 , $$ \begin{aligned} \begin{array}{cc} \dfrac{\dot{M}(\Omega )}{\dot{M}(\Omega = 0)} \approx \dfrac{(1-\Gamma )^{\frac{1}{\alpha }-1}}{\Big (1-\frac{4}{9} (\frac{v}{v_\mathrm{crit, 1} })^2 -\Gamma \Big )^{\frac{1}{\alpha }-1}}, \end{array} \end{aligned} $$(1)

where α is the force multiplier parameter, the critical velocity v crit , 1 = ( 2 3 GM R pb ) 0.5 $ v_{\mathrm{crit,1}} = \Big(\frac{2}{3}\frac{GM}{R_{\mathrm{pb}}}\Big)^{0.5} $, and Rpb is the polar radius at breakup, where we assume that the polar radius Rp does not alter with rotation rate (Ω), Rpb/Rp(Ω) = 1.

3. Nucleosynthesis

We present net wind yields and ejected masses for each Z grid of stellar models. Since stellar models with Z ≥ 0.1 Z⊙ include rotation, we can study the effects of rotation on the ejecta at these metallicities. Thus we compile a table of outputs for each Z, with the ZSMC and 0.1 Z⊙ grids including three subgrids with different rotation rates (0, 40% and 70% critically rotating). In line with Hirschi et al. (2005) and Higgins et al. (2023), the net wind yield calculated for a star of initial mass, m, and isotope, i, is

m i wind = ∫ 0 τ ( m ) M ˙ ( m , t ) [ X i S ( m , t ) − X i 0 ] d t , $$ \begin{aligned} m^\mathrm{wind}_{i} = \int _{0}^{\tau (m)} \dot{M}(m, t) \,[{X}^{S}_{i}(m, t) - {X}^0_{i}] \,dt, \end{aligned} $$(2)

where Ṁ is the mass-loss rate, X i S $ {X}^{S}_{i} $ is the surface abundance of a given isotope, and Xi0 is the initial abundance of a given isotope at the ZAMS. Net yields were then integrated from the beginning of core H burning until τ(m), corresponding to the end of core He burning for ZSMC and 0.1 Z⊙ grids and the end of core H burning for 1/30th Z⊙ and 0.01 Z⊙ grids, respectively.

We also calculated the ejected masses, EM, of each isotope, i, with

E M im = ∫ 0 τ ( m ) M ˙ X i S ( m , t ) d t . $$ \begin{aligned} EM_{im} = \int _{0}^{\tau (m)} \dot{M} \, {X}^{S}_{i}(m, t) \,dt . \end{aligned} $$(3)

The key reactions which dominate the wind yields of massive stars involve the CNO cycle, and secondary cycles (including the Ne-Na and Mg-Al cycles) during H burning which affect abundant isotopes of C, N, O, F, Ne, Na, Mg and Al. During core He burning, the H-synthesised 4He produces 12C via the triple-α reaction, resulting in an increase in 12C which activates the 12C(α, γ)16O reaction. The resulting CO core at core He exhaustion can subsequently be used to estimate the explodability of the stellar core (O’Connor & Ott 2011; Farmer et al. 2019).

At lower Z, the evolution of massive stars is hotter and more compact, relative to higher Z environments. Furthermore, the amount of mass lost during the core H burning phase decreases as a result of the Z-dependence on radiation-driven winds (see Fig. 1). In fact, the transition point at which massive and very massive stars switch from an optically thin to optically thick wind also shifts to higher initial masses at lower Z (Sabhahit et al. 2023). This means that ∼80–200 M⊙ stars will have optically thin winds at Z ∼ 0.1–0.01 Z⊙, and could result in more He-burning products in the ejecta rather than H-rich ejecta.

We explore the amount of mass lost in the H-burning phase compared to the He-burning phase, and what remnant mass is left for a range of initial masses at ZSMC and 0.1 Z⊙. As expected, the higher Z models with higher initial masses lose a significant proportion of their total mass already on the MS, with very little mass lost during the post-MS. Similarly, the remnant masses of the 0.1 Z⊙ models are significantly higher than ZSMC models with the same initial mass, because (i) the transition point has shifted from ∼100–200 M⊙ to ∼200–300 M⊙, and (ii) stars in lower Z environments lose less mass at all evolutionary stages. Interestingly, the same behaviour is seen from the He-burning ΔM at each Z and rotation rate, where lower mass stars lose more mass on the post-MS. With increasing initial mass, the post-MS winds become less relevant in comparison to both the MS winds and the remnant masses.

A key question bridging stellar winds and galactic chemical evolution is whether stellar yields are affected mainly by the Z-dependence on stellar winds or whether the nucleosynthesis is affected significantly enough to alter the ejecta. Figure 2 shows the surface evolution of isotopes (right to left, white region) and interior composition (grey shaded region) at core He exhaustion for a 100 M⊙ star at Z⊙, ZSMC, 10% Z⊙ and 1% Z⊙. Firstly, the dramatic change in the stellar mass at core He exhaustion is seen (grey region, marked by black line), also emphasising the increased ejected mass with increasing Z (white region). For instance, while the Z⊙ model ejects ∼20 M⊙ more than the 10% Z⊙ model, the 10% Z⊙ model loses approximately ten times more mass than the 1% Z⊙ model.

thumbnail Fig. 2.

Overview of the chemical composition and its evolution for a 100 M⊙ model calculated from the onset of core H burning until core He exhaustion for Z⊙ (upper left), Z = 0.004 (upper right), Z = 0.002 (lower left), and Z = 0.0002 (lower right). The shaded region illustrates the interior composition at the end of core He burning, while the unshaded region shows the surface abundances of each isotope with time evolving from right to left as mass is lost through winds. The vertical black solid line denotes the end of core He burning, while the vertical magenta solid line shows core H exhaustion.

Secondly, while the most abundant elements (H, He, C) are not altered in central mass fraction by Z, the less abundant elements such as Ne, Na, Mg and Al are shown to have different behaviours at lower Z. For instance, the amount of N lost on the MS at Z⊙ is higher than at ZSMC, leaving less N to be synthesised into 22Ne in the post-MS. One can see the change in 22Ne (orange dashed line) from the top panel to the lower panels. Of course, with decreased ZTotal each element will have a slightly lower initial abundance for each low Z environment, and different evolutionary timescales for each burning phase. However, by comparing the relative change in each isotope in the core (relative to their initial abundances) during the core H-burning phase as a function of increasing core He abundance (or MS evolution) we find that the abundances can change by ∼0.3 dex for 12C and 16O, and up to ∼ 1.5 dex for 23Na and 27Al, when comparing Z⊙ models to low Z models (Fig. 3).

thumbnail Fig. 3.

Central abundances of 12C (upper left), 16O (lower left), 23Na (upper right), and 27Al (lower right) relative to their initial abundance for each Z as a function of core He abundance during the MS evolution of 100 M⊙ models.

Finally, we examined the central abundance of 12C, 16O, 23Na and 27Al relative to their initial composition throughout the MS evolution (Yc) in Fig. 3. We find that 12C and 16O show central abundance patterns which are not significantly affected by Z during the MS, with only ∼0.2 dex difference between Z⊙ and 1% Z⊙ models. The central 12C quickly drops as a result of the CN-cycle, before slightly increasing again during the CNO-cycle and remains constant before a slight upturn towards the end of the MS. The Z⊙ model is most efficient at converting C to N, having the lowest relative abundance for the majority of the MS, but models of decreasing Z have higher C due to their higher central temperatures, which enable more rapid CNO burning. Similarly, the central 16O abundance drops at the expense of N at the onset of the MS, with the 1% Z⊙ model having the lowest central O abundance during the MS evolution. This confirms that the lower Z models have more efficient reactions during the MS due to their higher central temperatures. While the general behaviour of the more abundant CNO isotopes is similar for all Z, this is not always the case for heavier isotopes such as 23Na and 27Al. We find that all isotopes up to 23Na have the same general behaviour, yet from 23Na – 24, 25Mg – 26, 27Al there are differences in the reactions, mainly between Z > 10% Z⊙ and 1% Z⊙ models. For instance, the central 23Na abundance appears to increase steadily for the majority of the MS, but decreases at the end of the MS for ZSMC and 10% Z⊙ models.

Interestingly, we find that the central 23Na abundance drops dramatically in the 1% Z⊙ model, already in the first 10% of core H burning. Proton capture on 23Na proceeds through resonances at the higher core temperatures for 1% Z⊙ (Tc > 55 MK, see Fig. A.1), rapidly destroying the 23Na (Boeltzig et al. 2019). Moreover, the ratio between the 23Na(p,α)20Ne and the 23Na(p,γ)24Mg reaction changes from over a factor 100 to closer to a factor 10 at these temperatures enhancing the feeding to Mg and Al. Similarly, the central 27Al abundance increases slightly during the MS for Z > 1% Z⊙ models, but at 1% Z⊙ the behaviour is different with a dip mid-way through core H burning before steeply increasing until core H exhaustion.

The reason for these differing behaviours is a combination of the production of each isotope relative to the initial Z composition, which includes the abundance of the precursor isotopes in the Ne-Na and Mg-Al chains. The central temperatures are also slightly different at each Z, leading to changes in the efficiency of nuclear reactions. This means that the overall stellar properties may also play a role in the reaction rates and production of specific isotopes and net wind yields at lower Z. In general, the nuclear energy generation is more efficient at lower Z due to the increased central temperature and density, which has the largest effect on reactions that are heavily T-dependent. The luminosity generated by the CNO cycle, log10LCNO, attributed to the H-burning energy generation is increased from 6.2 L⊙ at Z⊙ to 6.3 at 1% Z⊙, though the largest difference in LCNO occurs in the remaining 20% of H burning where we see a 0.3 dex difference at Z⊙, which correlates with the steep increase in 27Al from Yc = 0.8–1.0.

While there are some changes to the nucleosynthesis of heavier elements at the lowest Z (1% Z⊙), the Z-dependence on winds may be the dominant effect on yields from massive stars. For instance, while the central abundances may change by 0.3 dex across Z (Z⊙–1% Z⊙) for the more abundant CNO isotopes, and up to 1 dex for the less abundant Ne-Na and Mg-Al chains, the mass-loss rates are 1 dex lower at 1% Z⊙ compared to Z⊙ during the first half of the MS, and up to 2 dex lower than Z⊙ models during the second half of the MS (see Fig. 1). Moreover, even if the central abundances changed significantly during the evolution of a low Z (1% Z⊙) star, if it is not lost in the winds before becoming reprocessed into another isotope then it will not affect the yields (see Fig. 2 bottom panel for example).

4. Rotation effects at low Z

In this section we present the wind yields of our ZSMC and 10% Z⊙ models which are calculated for three rotation rates (non-rotating, 40% critical and 70% critical). As discussed in Sect. 2, we implement the increase in mass-loss rates due to the proximity to the ΩΓ-limit as described by Maeder & Meynet (2000).

4.1. ZSMC

At 20% of Z⊙, the SMC is an excellent nearby laboratory for testing the evolution of massive stars at low Z. While the effects of stellar winds are reduced, VMSs (M > 100 M⊙) still lose 20–80% of their total mass on the MS and massive stars (30 M⊙ < M < 100 M⊙) lose 20–60% of their mass in the post-MS, leaving ∼10–40% of their initial mass. The Zenodo published tables provide the full list of ejected masses and net yields for isotopes of 1H to 28Si, from the ZAMS until core He exhaustion for three sets of critical rotation rates (0 top, 40% middle, and 70% bottom). We find that the amount of C, N, O, Ne, Na, and Mg ejected are still significant for all rotation rates (100 to 10−2 M⊙). We explore 100 M⊙ models at ZSMC for a range of rotation rates in Fig. 4. The Hertzsprung-Russell diagram (HRD) of these models (upper right) show that all 100 M⊙ models evolve bluewards after core H exhaustion at log Teff ∼ 4.0 and continue towards stripped helium stars forming cWRs. Interestingly, the mass-loss rates of the 40% critically rotating model are not sufficiently high enough to prevent spin up at the contraction phase, while the 70% critically rotating model loses enough mass and angular momentum to spin down throughout its evolution (upper left in Fig. 4). The surface enrichment of 14N increases with rotation, partly due to increased mass loss yet also due to increased rotational mixing dredging more 14N to the surface early on the MS. This illustrates that with increased rotation, stars will eject significantly more 14N throughout their MS evolution.

thumbnail Fig. 4.

Surface rotation (top left), HRD (top right), surface 14N evolution (bottom left) and mass evolution (bottom right) of 100 M⊙ models at ZSMC for a range of critical rotation rates (0% in solid red lines, 40% in dash-dotted green lines, and 70% in dotted grey lines).

Figure 5 shows the ejected mass of 14N from 30–500 M⊙ stars over the core H-burning phase for each rotation rate. Regardless of initial rotation rate (including non-rotating), the 200–500 M⊙ stars can contribute a significant amount of 14N in the range 0.1–1 M⊙ per individual VMS (Vink 2023). Below the transition point (M < 100 M⊙), where stars drive optically thin winds rather than the optically thick winds experienced by VMSs, the ejecta of 14N quickly drops. However, with rapid rotation, these 50–100 M⊙ stars can still eject a relevant quantity of 14N. When weighted by an IMF (M−1.9) the contribution from 200 M⊙ stars remains dominant (Fig. 6). Moreover, canonical O stars with masses 30–100 M⊙ require extreme rotation rates of Ωini/Ωcrit = 0.7 to compare with the ejecta of 300–500 M⊙ stars.

thumbnail Fig. 5.

Ejected mass of 14N as a function of initial mass, presented in solar masses, calculated at ZSMC during the core H-burning phase. Each coloured line illustrates a different rotation rate. Non-rotating models are shown in blue, 40% critically rotating models are in red, and 70% critically rotating models are in green.

thumbnail Fig. 6.

Initial mass function-weighted ejected masses of 14N as a function of initial mass, presented in solar masses, calculated at ZSMC for the core H-burning phase. Each coloured line illustrates a different rotation rate. Non-rotating models are shown in blue, 40% critically rotating models are in red, and 70% critically rotating models are in green.

When the total core H- and He-burning phases are considered, it is the rotating stars just below the transition point (50–80 M⊙) which dominate the IMF-weighted ejecta (Fig. 7). The total contribution from this mass range outweighs that of stars beyond 100 M⊙, as well as canonical O stars below 50 M⊙. Moreover, we show the IMF-weighted yields of 14N in Fig. 7 highlighting a peak in 14N ejecta from stars with initial masses 80–100 M⊙. The maximum 14N IMF-weighted net yields come from this mass range for all rotation rates, but particularly models with 40% critical rotation, while in the mass range of 100–300 M⊙, the 70% critically rotating models eject the most 14N for these masses.

thumbnail Fig. 7.

Initial mass function-weighted net wind yields of 14N as a function of initial mass, presented in solar masses, calculated for the H- and He-burning phases, at ZSMC. Each coloured line illustrates a different rotation rate. Non-rotating models are shown in green, 40% critically rotating in grey and 70% critically rotating models in orange.

In general, rotating massive stars lose more mass than non-rotating stars due to the proximity of the ΩΓ limit, which increases mass-loss rates due to the von Zeipel effect. Interestingly, at 70% critical rotation, we find that VMSs eject substantial amounts of C and O, ∼10 M⊙ more compared to the non-rotating and 40% critically rotating models. The non-rotating and 70% critically rotating models eject less 14N than the 40% rotating models since the rotational mixing effects lead to higher 14N enrichment at the surface earlier in the evolution. But the highest rotation rates also lead to higher mass-loss rates on the post-MS, thus ejecting more C and O. At the highest masses (∼200–500 M⊙), the rotation rate does not appear to have an effect on the yields of Ne, Na, Mg, Al, or Si. However, at the lower mass range ∼30–80 M⊙, rotating stars provide an order of magnitude more 14N than non-rotating stars when weighted by an IMF over the MS phase (Fig. 6).

4.2. 1/10th Z⊙

Figure 8 shows the stellar mass properties at both the end of core H burning and the end of core He burning as a function of initial stellar mass for models with two different metallicities: ZSMC (green lines) and 10% Z⊙ (black lines). The different line styles represent variations in the rotation rate, where solid lines correspond to zero rotation, dashed lines represent 40% of the critical rotation, and dash-dotted lines indicate 70% of the critical velocity. Across all panels, the most noticeable effect is that mass loss increases significantly with higher initial mass and metallicity. Stars at 10% Z⊙ lose less mass compared to those with ZSMC especially for stars with initial masses beyond ∼100 M⊙. Furthermore, the rotation rate significantly influences mass loss, with higher rotation rates (dash-dotted lines) leading to much lower final masses, especially at masses M < 300 M⊙. The indirect effects of rotation become important for stars with initial masses above ∼100 M⊙, where models with 70% of the critical rotation rate lose more mass due to the ΩΓ-limit early in the evolution. Subsequently stripping the star of its angular momentum, the resulting mass at core H exhaustion or final masses are much higher than stars which began their evolution with 40% critical rotation. The bottom left panel, which shows the same behaviour as the upper two panels demonstrate that rotation plays a role throughout the evolution by directly increasing the mass-loss rate, but also indirectly by increasing the luminosity and temperature of rotating models leading to increased mass loss even on the post-MS where most of the angular momentum has already been lost. By comparing the masses at core H exhaustion with final masses (top left and right panels) we see that stars still lose a significant amount of mass on the post-MS (see the gradient of the 1:1 mass ratio shown by grey dotted lines). The bottom right panel indicates that despite rotationally induced mass loss, the total mass lost as a function of initial mass remains relatively linear, though the 10% Z⊙ models lose up to 100 M⊙ less than ZSMC models at the highest mass range (∼400–500 M⊙).

thumbnail Fig. 8.

Stellar mass properties at the end of core H burning (MTAMS) and the end of core He burning (final) for models calculated at ZSMC (green) and 10% Z⊙ (black) for a range of critical rotation rates (0% in solid lines, 40% in dashed lines, and 70% in dash-dotted lines). A 1:1 mass ratio has been plotted for each panel in grey dotted lines.

Figures 9 and 10 show the effects of rotation on the evolution of 50 M⊙ and 300 M⊙ stars at 10% Z⊙. The 50 M⊙ models (Fig. 9) follow a standard O-star mass loss prescription, while the 300 M⊙ models (Fig. 10) are subject to enhanced optically thick winds, resulting in notable differences in mass loss and rotational evolution.

thumbnail Fig. 9.

Surface rotation (top left), HRD (top right), surface 14N evolution (bottom left), and mass evolution (bottom right) of 50 M⊙ models at 10% Z⊙ for a range of critical rotation rates (0% in solid red lines, 40% in dash-dotted green lines, and 70% in dotted grey lines).

thumbnail Fig. 10.

Surface rotation (top left), HRD (top right), surface 14N evolution (bottom left), and mass evolution (bottom right) of 300 M⊙ models at 10% Z⊙ for a range of critical rotation rates (0% in solid red lines, 40% in dash-dotted green lines, and 70% in dotted grey lines).

In the 50 M⊙ models, rotation plays a moderate role in altering surface properties. As shown in the upper left panel, higher rotation rates (Ωini/Ωcrit = 40% and 70%) maintain relatively high surface velocities for the entire H-burning phase. By the end of core He burning, even the most rapidly rotating models exhibit significant angular momentum loss as a result of wind mass loss. The HRD (upper right panel) reveals that the faster rotating models shift to slightly higher luminosities and hotter effective temperatures as a result of rotational mixing. However, the O-star wind rates limit the overall loss of angular momentum, leading to significant surface 14N enrichment, and relatively minor differences in mass evolution (for non-rotating and 40% critical models). The 70% critical rotation model experiences significant chemical mixing and evolves bluewards to hotter temperatures towards the end of core He burning and loses a substantial amount of mass in this burning phase.

In contrast, the 300 M⊙ models display a much more significant impact of rotationally induced mass loss since these stars are above the transition point and already experience enhanced optically thick winds with high mass-loss rates. Hence, the combined effect of rotation and mass loss can play a key role in the evolution of stars near the transition point at low metallicity. As seen in the surface velocity evolution (upper left panel), stars with higher rotation rates experience a rapid decline in surface velocity from the onset of core H burning, which is driven by the strong mass loss from optically thick winds. This results in the 70% rotating model losing a considerable portion of its angular momentum early on, leading to a rapid decrease in surface rotation. The HRD (upper right panel) shows that rotating models evolve closer to the transition point causing the switch from horizontal evolution (non-rotating) to the vertical evolution of VMSs seen at predominantly higher Z, with a dramatic drop in luminosity and subsequent WR evolution. The mass evolution plot (bottom right panel) highlights the impact of rotationally enhanced winds, with the 70% rotating model experiencing significant mass loss during the MS, leading to a 50–100 M⊙ difference with the non-rotating model during the MS. Additionally, the surface 14N enrichment is more pronounced in the 300 M⊙ model (bottom left panel), with rapidly rotating models showing a very quick increase by a factor of 10 in 14N abundance. This rapid increase in 14N is due to higher central temperatures (than lower mass or higher Z stars) reaching CN-equilibrium almost instantly, stripping of the outer pristine layers quickly, and the effects of rotational mixing dredging 14N to the surface. Comparable Figures A.2 and A.3 for 50 M⊙ and 300 M⊙ models at ZSMC are shown in Appendix A. These figures illustrate that at higher Z, the effects of rotation are diminished slightly given the increased mass loss, particularly for the 300 M⊙ models which are now far beyond the transition point.

5. Globular clusters

In this section we address the importance of GCs at Z ∼ 3% Z⊙ for stellar nucleosynthesis and evolution. Such low metallicity GCs are common in the outer regions of galaxies, including our Milky Way, and are considered to be remnants from the early universe, providing a glimpse into the conditions of early star formation. One of the most significant chemical peculiarities observed in GCs is the Na-O anti-correlation. This refers to the inverse relationship between the abundance of sodium (Na) and oxygen (O) in stars within the same cluster. Stars with high Na content typically have low O content, and vice versa. This pattern is not seen in field stars, making it unique to GCs, and is expected to be caused by multiple stellar populations. These populations are characterised by different chemical compositions, including the Na-O anti-correlation. Population 1 (P1) stars generally show abundances similar to field stars of similar metallicity, while Population 2 (P2) stars exhibit enhanced Na and depleted O, suggesting a different chemical enrichment history. The distinction between P1 and P2 stars suggests a complex formation history for GCs. It implies that clusters formed multiple populations of stars over time, with the second population forming from material that was chemically enriched.

Carretta et al. (2009) explore the Na-O anti-correlation in 15 GCs and the distinction between different stellar populations within these clusters. They define P1 stars as those with similar [Na/Fe] and [O/Fe] ratios as field stars, while stars with enhanced Na/Fe and depleted O/Fe are considered P2 stars. We investigate the relevance of our 80–500 M⊙ models in forming P2 stars by exploring the surface evolution of [Na/Fe] and [O/Fe] in relation to the observations of stars in NGC 5904. While this cluster is not atypical for GCs, it provides a broad range of [Na/Fe] and [O/Fe] to test our evolutionary models against. Moreover, the UVES observed [Fe/H] abundance of NGC 5904 (−1.34) is close to our model grid (−1.35). We note that since our initial composition does not account for the α-enhancement observed (e.g. Gratton et al. 2012), we have corrected our model data in [O/Fe] by increasing the ratio at each point by +0.4 dex. Figure 11 illustrates the MS evolution of the surface composition in our models with respect to the abundance ratios of observed stars in NGC 5904. We find that stars with Mi ≥ 200 M⊙ qualitatively explain the P2 stars lying outside the green circle which shows the abundances of field stars (P1) where [O/Na] ≃ 0.6 (Charbonnel 2016). As the models evolve during core H burning, their surfaces enrich further in Na and deplete in O. All models begin with the same abundance ratios relative to solar (0, 0). The higher mass models become more O-depleted than lower mass models, and also display increased Na-enrichment.

thumbnail Fig. 11.

Surface evolution of [Na/Fe] and [O/Fe] for 80–500 M⊙ models during core H burning at 1/30th Z⊙. The red circles represent detections of stars in NGC 5904, while blue triangles represent upper limits in [O/Fe] as defined by Carretta et al. (2009). The green dashed circle establishes the abundance ratios of field stars ([O/Na] = 0.6), with stars outside of it representing P2 stars.

Similarly, we investigate the surface evolution of Na-O during the post-MS in Fig. 12 finding that stars with Mi ∼ 80 − 100 M⊙ evolve with surface Na and O ratios that encompass the observed Na/O pattern. These figures illustrate the effects of winds on the surface evolution and ejecta of M > 80 M⊙ stars. As a result of stronger mass loss on the MS, higher mass models (M > 200 M⊙) lie beyond the observed [Na/O] pattern, with extremely O-depleted surfaces until both O and Na enrich as a result of He-burning nucleosynthesis. However, even when these VMS models enrich in O during core He burning, the relative [O/Fe] abundances are orders of magnitude higher than observed P2 stars. Therefore, the best representation of observed NGC 5904 abundances appears to be VMSs with M ≥ 200 M⊙ on the MS or 80–100 M⊙ stars on the post-MS. At low Z (3% Z⊙), MS winds are reduced for stars below the transition point (M ≤ 200 M⊙) and therefore do not showcase the H-processed Na and O until much later in their evolution. On the other hand, VMSs (M ≥ 200 M⊙) strip their outer layers more quickly to reveal the Na-rich and O-depleted material early on the MS.

thumbnail Fig. 12.

Ratio of [Na/Fe] as a function of [O/Fe] for 80–500 M⊙ models during core He burning at 1/30th Z⊙. First and second population stars from the cluster NGC5904, adopted from Carretta et al. (2009), are shown in red circles and blue triangles, respectively.

Figures 11 and 12 illustrate the evolution of surface abundances of our models during the MS and post-MS, which provides a comparison to the observed abundance ratios in GCs such as NGC 5904. While these surface properties highlight the relevance of VMSs in reproducing the abundance patterns of Na and O in GCs, they do not directly correlate to the ejecta of these stars which could form a second population of stars. Therefore, we also present the net wind yields of 16O and 23Na, relative to their initial composition, for each of our models integrated over the H-burning phase. Table 2 confirms that on the MS, stars with Mi ≥ 100 M⊙ eject O-depleted, Na-rich material with up to 10−3 M⊙ of Na ejected for the highest mass models.

Table 2.

Net wind yields of 16O and 23Na for models calculated at 3% Z⊙ and initial masses of 30–500 M⊙ in solar masses calculated from the onset of core H burning until core H exhaustion.

This scenario changes for VMSs when the post-MS is included in the total net yields of Na and O, integrating the abundances of mass lost from the onset of core H burning until core He exhaustion. At this point, sufficient O-rich material is ejected (10−4 M⊙) to result in positive net yields of O. Therefore, we find that the MS yields of VMSs best represent the observed trend in Na-O abundances. Crucially, at such low Z, the net wind yields of stars with M < 100 M⊙ are 10 orders of magnitude lower than stars above this mass range, suggesting that massive O stars cannot produce enough Na-rich and O-poor material compared to VMSs to produce a second population of stars with such abundance ratios.

Another way to consider the contribution of [Na/Fe] enhanced and [O/Fe] depleted material from these stellar models is to investigate the interior and surface composition in mass fraction with time. Figures 13 and 14 show the same surface composition evolution as in Fig. 2, where the interior composition remaining at the end of core He burning is shown in the grey shaded region. Though at this point, we also include the [Na/Fe] and [O/Fe] ratios in pink and green respectively, starting at 0 (right) and evolving (left) as mass is lost from the outer layers. One can see in Fig. 13 that there is a small region (∼72 M⊙) where the surface is [Na/Fe] rich, but it is only after the black solid line, which represents the end of core He burning. Within the grey shaded interior composition from surface to centre (right to left), there is a significant drop in [O/Fe] with a simultaneous rise in [Na/Fe] which then remains constant in the region 35–60 M⊙ in the stellar interior. This means that during 99% of the stellar lifetime, an 80 M⊙ star does not eject a large amount of material (at this Z) which matches the observed Na-O pattern, but if the star did experience enhanced mass loss towards the end of its life through eruptions or another mechanism, then the interior composition which is representative of the Na-O pattern could be lost to the ISM. On the other hand, a 400 M⊙ star with the same Z and initial composition shown in Fig. 14 produces a similar portion of [Na/Fe]-rich and [O/Fe]-poor material (see region ∼220–370 M⊙). However, the plateau in both abundance ratios occurs over a larger mass range (∼150 M⊙ compared to ∼30 M⊙ in the 80 M⊙ model) and crucially, is lost in the ejecta (white region) rather than being maintained in the stellar interior (grey region).

thumbnail Fig. 13.

Time evolution of surface abundances of a 80 M⊙ model at 3% Z⊙ calculated from the onset of core H burning until core He exhaustion. The [Na/Fe] and [O/Fe] are included.

thumbnail Fig. 14.

Time evolution of surface abundances of a 400 M⊙ model at 3% Z⊙ calculated from the onset of core H burning until core He exhaustion. The [Na/Fe] and [O/Fe] are included.

Finally, we showcase the ejecta from our models alongside observed abundance patterns to directly compare the composition of P2 stars with the composition of the ejected mass from 80–500 M⊙ stars with the same [Fe/H] metallicity. Figure 15 presents the cumulative ejected masses integrated from the onset of core H burning until core He exhaustion by adopting equation (3) and integrating over short timescales covered in ten model steps throughout the entire evolution. After calculating the ejected mass of Na and O for each segment, the ejecta are normalised to the [Na/Fe] and [O/Fe] notation, by scaling the log (Na/Fe) fractions by number relative to solar values. Each integrated ejected mass is presented by a marker, which is different for each initial mass model, and then connected by a dashed line to show the evolution of these cumulative ejected masses over time as they become increasingly Na rich and O depleted. The same observations from NGC 5904 shown in Figs. 11 and 12 are included, where the green ellipse encompasses P1 stars which have a surface composition representative of field stars, while those lying beyond the green ellipse are considered P2 stars. Figure 15 confirms that all models begin their evolution by ejecting material with a composition (relative to solar) similar to P1 stars. Then, as the higher mass models begin to peel their outer layers to expose fresh nucleosynthesised material, the ejecta quickly becomes Na rich and O poor. This behaviour occurs for all masses in the range 80–500 M⊙ but as previously discussed, happens later in the evolution for 80–100 M⊙ models than for 200–500 M⊙ models. Therefore, the longest evolutionary phase with the highest mass loss, which provides ejected masses with abundance patterns which are representative of the P2 stars, is the MS evolution of VMSs with masses ∼200 M⊙.

thumbnail Fig. 15.

Cumulative ejected masses of [Na/Fe] and [O/Fe] for 80–500 M⊙ models calculated at 1/30th Z⊙. Observations of stars in NGC 5904 are shown as in Fig. 11, adopted from Carretta et al. (2009). The insert includes the wide range of O-enrichment ejected by the 400–500 M⊙ models during the post-MS.

A key point in establishing the second population in GCs is that the ejected material should be slow-moving such that it remains in the potential well of the cluster to form a second population of stars. Supernovae ejecta and WR winds experience velocities of a few thousand km/s which are much higher than the predicted escape velocities of GCs (e.g. Gieles et al. 2018), while VMS winds are likely slower than those from WR stars (Vink 2018, 2023). In fact, since VMSs can eject synthesised material from the onset of core H burning for a few million years, and the second population of stars in GCs can be formed within a few million years of the first population (Gratton et al. 2012), VMSs could in fact be an important source for forming P2 stars in GCs (Lochhaas & Thompson 2017). Moreover, Gratton et al. (2012) discusses the observed spread in Fe abundances in Galactic GCs which resembles similar patterns observed in local dwarf galaxies and suggests that GCs may start their evolution as the cores of dwarf galaxies (Bekki & Freeman 2003; Bekki et al. 2007; Böker 2008; Carretta et al. 2010). Interestingly, a large homogeneous sample of over 900 massive stars in the 30 Doradus region of the Large Magellanic Cloud was investigated by the VLT-FLAMES Tarantula Survey (VFTS; Evans et al. 2011), exploring the complex star-formation of the nebula, as well as detecting the most massive stars to date with M* ∼ 200 M⊙. Furthermore, Schneider et al. (2018) found an excess of massive stars in the region and a top-heavy IMF with a power law index of 1.90. These recent studies collectively adhere to the theory that VMSs may be the most efficient solution to self-enrichment in GCs (Vink 2018).

6. Discussion

6.1. Dilution

The extent of the Na-O anti-correlation observed across numerous clusters has provided important clues for constraining the multiple population phenomenon in GCs. While some ejecta from the first population may be encapsulated to form the second population, it is likely that a fraction of the mass accumulated will be from the original cluster gas. Therefore, we tested the effect of diluting our wind ejecta with the initial composition to estimate the dilution factor required to reproduce the observed abundance ratios of Na and O. We calculated the dilution as

X D = X i × f − 1 f + M EM Δ M × 1 f , $$ \begin{aligned} X_{D} = X_{i} \times \frac{f-1}{f} + \frac{M_{EM}}{\Delta M} \times \frac{1}{f}, \end{aligned} $$(4)

where XD is the abundance of the diluted ejecta, Xi is the initial abundance of the original gas for a given isotope, f is the dilution factor, MEM is the total wind ejecta of a given isotope, and ΔM represents the total mass lost. Figure 16 shows the total cumulative ejecta of [Na/Fe] and [O/Fe] for models with initial masses ranging from 80–300 M⊙ since stars below 80 M⊙ do not eject a meaningful amount of mass, and stars beyond 300 M⊙ eject too much O-rich material to provide a total ejecta representative of the O-depleted population. We dilute the ejecta by factors ranging from 1 to 100, and show the resulting diluted gas abundance for each f value with a different coloured line. We find that we require a low dilution factor of f = 2–4 in order to best represent the observed [Na/Fe] abundance pattern, since VMS with M = 100 − 300 M⊙ can reach the increased [Na/Fe] values of up to 0.8 dex. The [O/Fe] ratio is more challenging to explain with diluted values. However, the undiluted VMS models can reach the negative [O/Fe] values, thus possibly explaining the full range of observed abundances.

thumbnail Fig. 16.

Diluted ejected from 80–500 M⊙ models calculated at 1/30th Z⊙. Observations of stars in NGC 5904 are shown as in Fig. 11, adopted from Carretta et al. (2009). Dilution factors ranging from 1 to 100 are shown for each model (assorted symbols) and connected for each f value (coloured lines).

We included a subset of rotating models calculated at 1/30th Z⊙, with Ωini/Ωcrit = 0.4, and we note that some models do not complete the core He burning since higher rotation rates lead to stars approaching criticality at this low Z. The impact of rotation on the observed [Na/Fe]-[O/Fe] correlation is shown in Fig. 17 for a range of dilution factors, comparable to Fig. 16. We find that the 80 M⊙ model does not lose significant mass on the main sequence but loses significant amounts of 16O during core He-burning, which results in O-production rather than depletion regardless of dilution factor. The 80 M⊙ rotating model reaches the cWR phase at effective temperatures ranging from log Teff = 5.0–5.2, leading to 0.1 M⊙ of 16O ejecta. This is in contrast to the non-rotating models shown in Fig. 16 which displays a trend of increasing Na-enrichment and O-depletion. The 100–300 M⊙ models are not affected in the same way since their convective cores are almost the entire stellar mass, leading to negligible effects of rotational mixing and therefore minor changes in Na production and O destruction. Figure 17 illustrates that Na is overproduced in all rotating models, and by diluting the material, the [Na/Fe] shifts towards lower values until reaching the initial composition representative of field stars with f = 100. On the other hand, the difference between O depletion by the 100–300 M⊙ models and O production by the 80 M⊙ model leads to 2 paths of diluting back to the initial composition, one which ranges from ∼−1.0 up to 0.4 (similar to the behaviour seen in Fig. 16), or from 1.5 down to the original gas value of 0.4, respectively.

thumbnail Fig. 17.

Diluted ejected from initially 40% critically rotating 80–300 M⊙ models calculated at 1/30th Z⊙. Observations of stars in NGC 5904 are shown as in Fig. 11, adopted from Carretta et al. (2009). Dilution factors ranging from 1 to 100 are shown for each model (assorted symbols) and connected for each f value (coloured lines).

6.2. Mg-Al anti-correlation

The Na-O anti-correlation observed in GCs has been the focus of many studies (e.g. Carretta et al. 2009, 2010), with many progenitors explored to reproduce the observed abundance patterns and ultimately explain the formation of multiple populations. Another indicator of the second population in GCs is the observed Mg-Al anti-correlation, which has not been observed in all clusters, and in many cases is not as extended or evident in less massive or metal-rich clusters. In some sense, it is clear that when considering the Na-O anti-correlation, we expect Na to enrich with nucleosynthesis and O to be depleted (at least in the early stages of nuclear burning where most of the stellar lifetime is spent). In comparison, the Mg-Al cycle replenishes 24Mg while also forming Al, and therefore, it is less evident how such Mg-depleted stars could have such Al enrichment. In Fig. 18 we present our non-rotating model grid calculated at 1/30th Z⊙ alongside observations from NGC 7089 (M2) which has a similar [Fe/H] abundance (−1.47 dex ± 0.03) and a broad range of Mg and Al. We find that our models (scaled for α-enhancement) can represent the wide range of Al-enhancement observed, when considering the highest mass models. On the other hand, even when including all isotopes which contribute to observed Mg (24, 25, 26Mg and 26Al which decays to 26Mg), the models do not become sufficiently Mg-depleted to enclose the observed data. Decressin et al. (2007) suggested that the 24Mg(p, γ) reaction would need to be increased by a factor of 1000 to reproduce the observed Mg-depletion. Yet while many experimental rates suggest an upper/lower limit corresponding to ∼20% uncertainty in this rate, it is unlikely that the reaction would be revised to such an extreme case. Moreover, Choplin et al. (2018) suggests a similar increase by a factor 100 in the 23Na(p, γ)24Mg reaction, yet a significant change to this reaction rate would also effect the Na-enrichment observed in GCs due to the depletion of 23Na at the expense of 24Mg.

thumbnail Fig. 18.

Cumulative ejected masses of [Mg/Fe] and [Al/Fe] for 80–500 M⊙ models calculated at 1/30th Z⊙. Observations of stars in M2 are shown by blue triangles. The [Mg/Fe] ratio includes isotopes of 24, 25, 26Mg and 26Al, while the [Al/Fe] ratio is calculated for 27Al.

6.3. Wind yields at low Z

The contribution of nuclear-processed elements from stellar winds at low Z is considered to be diminished due to the effect of the metallicity dependence on radiatively driven winds. This may be the case for canonical O stars with initial masses ranging from 20–80 M⊙ which drive optically thin winds with mass-loss rates of order 10−6 to 10−8. However, rapid rotation can induce higher mass-loss rates as stars get closer to their ΩΓ-limit, leading to increased winds in lower Z environments. In this scenario, massive stars would need to rotate at rapid rotation rates (Ωini/Ωcrit > 0.4) in order to increase their wind rate sufficiently to impact their overall wind yields. We note that while the effect of rotational mixing at low metallicity is an interesting topic of stellar evolution, we do not investigate the full effects of rotation in this work. Since VMS are almost fully convective, with cores which are ≥90% of the total stellar mass, the impact of rotational mixing at these high masses are significantly reduced, even at the lowest Z. However, we tested the effect of rotationally enhanced winds due to the ΩΓ-limit as a function of chemical wind yields, but we note that at Z < 10% Z⊙, models with Ωini/Ωcrit > 0.4 reach criticality early in the evolution. Furthermore, what occurs in nature to resolve stars at criticality is undetermined, and potentially points to unresolved physics which is beyond the scope of this work.

On the other hand, above ∼80 M⊙ stars approach the transition point where more massive stars launch optically thick winds which lead to enhanced mass-loss rates of up to 10−3 even at low Z. This means that while the base mass-loss rates are reduced with Z, VMS above the transition point can still be impactful in their wind ejecta. We find that at 10% Z⊙ stars with initial masses ranging from 200–500 M⊙ still lose 60–90% of their total mass, and at 1/30th Z⊙ the most massive stars continue to lose up to 60% of their total mass, ejecting hundreds of solar masses of nucleosynthesised material into the ISM (see Zenodo tables). This suggests that VMS maintain a key role in the enrichment of the early Universe, particularly in their pollution of H-processed elements such as 14N and 23Na.

7. Conclusions

In this work, we have explored the effect of enhanced optically thick winds on the nucleosynthesis and chemical yields of massive and very massive stars across low metallicities. We provided a range of initial masses for four metallicities scaling from ZSMC down to 1% Z⊙. We tested the implementation of rotationally enhanced mass loss due to the ΩΓ limit for ZSMC and 10% Z⊙ models and focused on non-rotating models for 3% and 1% Z⊙ models. We provide ejected mass and net wind yields for all grids of models in Appendix B, which can be found at Zenodo and arXiV.

We find that due to the increasingly hot and compact evolution of stars with lower Z, the central temperature changes can have an effect on the nucleosynthesis of light elements from Na-Mg-Al isotopes. However, while the relative changes in these reactions and subsequent surface abundances may change by 0.1–0.3 dex, the mass-loss rates can change by orders of magnitude across this Z range. Therefore, the wind rates dominate the yields of a given isotope instead of the effect of the host Z environment on the nucleosynthesis.

Rotating (70% critical) models can lose up to 100 M⊙ more than their non-rotating counterparts during the MS. The effects of rotation are diminished at the highest mass range due to the loss of angular momentum early in the evolution. Counter-intuitively, the indirect effects of mass loss on rotation, while also considering the effects of rotationally induced mass loss, can mean that intermediate rotation rates (40% critical) can lose the most mass, as they spin up at low Z towards critical rotation, which results in higher ejecta masses than higher or lower rotation rates.

A Na-O anti-correlation is produced at the surface of VMSs (M > 100 M⊙) during the MS as well as at the surface of 80–100 M⊙ stars during the post-MS phase. Furthermore, the substantial amounts of Na-rich and O-depleted material lost by 200–300 M⊙ mass stars on their MS suggests that this phase would be most effective at enriching the ISM to form a second population of stars. We compared the normalised ejected masses of our 3% Z⊙ models with observations of red giant stars from the GC NGC 5904. Without tailoring model inputs to fit the observations, we find that our models naturally produce Na-enriched and O-depleted material, qualitatively reproducing the observed trend of [Na/Fe] and [O/Fe] with the surface abundances of the second population.

Data availability

Tables of ejected masses and stellar properties are available to download at: https://zenodo.org/records/15283909.

Acknowledgments

The authors acknowledge MESA authors and developers for their continued revisions and public accessibility of the code. JSV, AML, and ERH are supported by STFC funding under grant number ST/V000233/1 in the context of the BRIDGCE UK Network. RH acknowledges support from STFC, the World Premier International Research Centre Initiative (WPI Initiative), MEXT, Japan and the IReNA AccelNet Network of Networks (National Science Foundation, Grant No. OISE-1927130). This article is based upon work from the ChETEC COST Action (CA16117) and the European Union’s Horizon 2020 research and innovation programme (ChETEC-INFRA, Grant No. 101008324).

References

  1. Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197 [Google Scholar]
  2. Arnett, W. D., Meakin, C., Hirschi, R., et al. 2019, ApJ, 882, 18 [NASA ADS] [CrossRef] [Google Scholar]
  3. Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
  4. Baraffe, I., Heger, A., & Woosley, S. E. 2001, ApJ, 550, 890 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bastian, N., & Lardo, C. 2018, ARA&A, 56, 83 [Google Scholar]
  6. Bekki, K., & Freeman, K. C. 2003, MNRAS, 346, L11 [Google Scholar]
  7. Bekki, K., Campbell, S. W., Lattanzio, J. C., & Norris, J. E. 2007, MNRAS, 377, 335 [Google Scholar]
  8. Boeltzig, A., Best, A., Pantaleo, F. R., et al. 2019, Phys. Lett. B, 795, 122 [CrossRef] [Google Scholar]
  9. Böker, T. 2008, ApJ, 672, L111 [CrossRef] [Google Scholar]
  10. Bowman, D. M. 2020, Front. Astron. Space Sci., 7, 70 [Google Scholar]
  11. Bromm, V., & Larson, R. B. 2004, ARA&A, 42, 79 [NASA ADS] [CrossRef] [Google Scholar]
  12. Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Carretta, E., Bragaglia, A., Gratton, R. G., et al. 2009, A&A, 505, 117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Carretta, E., Bragaglia, A., Gratton, R. G., et al. 2010, ApJ, 714, L7 [Google Scholar]
  15. Charbonnel, C. 2016, EAS Publ. Ser., 80-81, 177 [CrossRef] [EDP Sciences] [Google Scholar]
  16. Choplin, A., Meynet, G., Maeder, A., Hirschi, R., & Chiappini, C. 2018, J. Phys. Conf. Ser., 940, 012021 [Google Scholar]
  17. Cristini, A., Meakin, C., Hirschi, R., et al. 2017, MNRAS, 471, 279 [NASA ADS] [CrossRef] [Google Scholar]
  18. Cristini, A., Hirschi, R., Meakin, C., et al. 2019, MNRAS, 484, 4645 [NASA ADS] [CrossRef] [Google Scholar]
  19. Crowther, P. A., Caballero-Nieves, S. M., Bostroem, K. A., et al. 2016, MNRAS, 458, 624 [Google Scholar]
  20. Cyburt, R. H., Amthor, A. M., Ferguson, R., et al. 2010, ApJS, 189, 240 [NASA ADS] [CrossRef] [Google Scholar]
  21. Decressin, T., Meynet, G., Charbonnel, C., Prantzos, N., & Ekström, S. 2007, A&A, 464, 1029 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. de Jager, C., Nieuwenhuijzen, H., & van der Hucht, K. A. 1988, A&AS, 72, 259 [NASA ADS] [Google Scholar]
  23. Denissenkov, P. A., & Hartwick, F. D. A. 2014, MNRAS, 437, L21 [Google Scholar]
  24. Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146 [Google Scholar]
  25. Evans, C. J., Taylor, W. D., Hénault-Brunet, V., et al. 2011, A&A, 530, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Farmer, R., Renzo, M., de Mink, S. E., Marchant, P., & Justham, S. 2019, ApJ, 887, 53 [Google Scholar]
  27. Freytag, B., Ludwig, H. G., & Steffen, M. 1996, A&A, 313, 497 [NASA ADS] [Google Scholar]
  28. Gieles, M., Charbonnel, C., Krause, M. G. H., et al. 2018, MNRAS, 478, 2461 [Google Scholar]
  29. Gratton, R., Sneden, C., & Carretta, E. 2004, ARA&A, 42, 385 [Google Scholar]
  30. Gratton, R. G., Carretta, E., & Bragaglia, A. 2012, A&A Rev., 20, 50 [NASA ADS] [CrossRef] [Google Scholar]
  31. Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368 [NASA ADS] [CrossRef] [Google Scholar]
  32. Herwig, F. 2000, A&A, 360, 952 [NASA ADS] [Google Scholar]
  33. Higgins, E. R., & Vink, J. S. 2019, A&A, 622, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Higgins, E. R., Vink, J. S., Hirschi, R., Laird, A. M., & Sabhahit, G. N. 2023, MNRAS, 526, 534 [NASA ADS] [CrossRef] [Google Scholar]
  35. Hirschi, R., Meynet, G., & Maeder, A. 2005, A&A, 433, 1013 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Krause, M., Fierlinger, K., Diehl, R., et al. 2013, A&A, 550, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Lochhaas, C., & Thompson, T. A. 2017, MNRAS, 470, 977 [Google Scholar]
  38. Maeder, A., & Meynet, G. 2000, A&A, 361, 159 [NASA ADS] [Google Scholar]
  39. O’Connor, E., & Ott, C. D. 2011, ApJ, 730, 70 [Google Scholar]
  40. Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
  41. Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [Google Scholar]
  42. Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
  43. Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34 [NASA ADS] [CrossRef] [Google Scholar]
  44. Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10 [Google Scholar]
  45. Pols, O. R., Tout, C. A., Eggleton, P. P., & Han, Z. 1995, MNRAS, 274, 964 [Google Scholar]
  46. Rizzuti, F., Hirschi, R., Georgy, C., et al. 2022, MNRAS, 515, 4013 [NASA ADS] [CrossRef] [Google Scholar]
  47. Rogers, F., & Nayfonov, A. 2002, ApJ, 576, 1064 [NASA ADS] [CrossRef] [Google Scholar]
  48. Sabhahit, G. N., Vink, J. S., Sander, A. A. C., & Higgins, E. R. 2023, MNRAS, 524, 1529 [NASA ADS] [CrossRef] [Google Scholar]
  49. Sander, A. A. C., & Vink, J. S. 2020, MNRAS, 499, 873 [Google Scholar]
  50. Schaerer, D., Guibert, J., Marques-Chaves, R., & Martins, F. 2025, A&A, 693, A271 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Schneider, F. R. N., Sana, H., Evans, C. J., et al. 2018, Science, 359, 69 [NASA ADS] [CrossRef] [Google Scholar]
  52. Scott, L. J. A., Hirschi, R., Georgy, C., et al. 2021, MNRAS, 503, 4208 [NASA ADS] [CrossRef] [Google Scholar]
  53. Vink, J. S. 2018, A&A, 615, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Vink, J. S. 2023, A&A, 679, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Vink, J. S., & Harries, T. J. 2017, A&A, 603, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Vink, J. S., Muijres, L. E., Anthonisse, B., et al. 2011, A&A, 531, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  58. Vink, J. S., Heger, A., Krumholz, M. R., et al. 2015, Highlights Astron., 16, 51 [Google Scholar]
  59. Yusof, N., Hirschi, R., Meynet, G., et al. 2013, MNRAS, 433, 1114 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Additional figures

thumbnail Fig. A.1.

Central temperature evolution for 100 M⊙ models during core H burning calculated at various Z.

thumbnail Fig. A.2.

Surface rotation (top left), HRD (top right), surface 14N evolution (bottom left), and mass evolution (bottom right) of 50 M⊙ models at ZSMC for a range of critical rotation rates (0% in solid red lines, 40% in dash-dotted green lines, and 70% in dotted grey lines).

thumbnail Fig. A.3.

Surface rotation (top left), HRD (top right), surface 14N evolution (bottom left), and mass evolution (bottom right) of 300 M⊙ models at ZSMC for a range of critical rotation rates (0% in solid red lines, 40% in dash-dotted green lines, and 70% in dotted grey lines).

All Tables

Table 1.

Initial composition of our model grid for a range of initial Z in mass fractions.

In the text
Table 2.

Net wind yields of 16O and 23Na for models calculated at 3% Z⊙ and initial masses of 30–500 M⊙ in solar masses calculated from the onset of core H burning until core H exhaustion.

In the text

All Figures

thumbnail Fig. 1.

Mass-loss rates of 100 M⊙ models for each Z as a function of core He abundance during the MS evolution.
In the text
thumbnail Fig. 2.

Overview of the chemical composition and its evolution for a 100 M⊙ model calculated from the onset of core H burning until core He exhaustion for Z⊙ (upper left), Z = 0.004 (upper right), Z = 0.002 (lower left), and Z = 0.0002 (lower right). The shaded region illustrates the interior composition at the end of core He burning, while the unshaded region shows the surface abundances of each isotope with time evolving from right to left as mass is lost through winds. The vertical black solid line denotes the end of core He burning, while the vertical magenta solid line shows core H exhaustion.
In the text
thumbnail Fig. 3.

Central abundances of 12C (upper left), 16O (lower left), 23Na (upper right), and 27Al (lower right) relative to their initial abundance for each Z as a function of core He abundance during the MS evolution of 100 M⊙ models.
In the text
thumbnail Fig. 4.

Surface rotation (top left), HRD (top right), surface 14N evolution (bottom left) and mass evolution (bottom right) of 100 M⊙ models at ZSMC for a range of critical rotation rates (0% in solid red lines, 40% in dash-dotted green lines, and 70% in dotted grey lines).
In the text
thumbnail Fig. 5.

Ejected mass of 14N as a function of initial mass, presented in solar masses, calculated at ZSMC during the core H-burning phase. Each coloured line illustrates a different rotation rate. Non-rotating models are shown in blue, 40% critically rotating models are in red, and 70% critically rotating models are in green.
In the text
thumbnail Fig. 6.

Initial mass function-weighted ejected masses of 14N as a function of initial mass, presented in solar masses, calculated at ZSMC for the core H-burning phase. Each coloured line illustrates a different rotation rate. Non-rotating models are shown in blue, 40% critically rotating models are in red, and 70% critically rotating models are in green.
In the text
thumbnail Fig. 7.

Initial mass function-weighted net wind yields of 14N as a function of initial mass, presented in solar masses, calculated for the H- and He-burning phases, at ZSMC. Each coloured line illustrates a different rotation rate. Non-rotating models are shown in green, 40% critically rotating in grey and 70% critically rotating models in orange.
In the text
thumbnail Fig. 8.

Stellar mass properties at the end of core H burning (MTAMS) and the end of core He burning (final) for models calculated at ZSMC (green) and 10% Z⊙ (black) for a range of critical rotation rates (0% in solid lines, 40% in dashed lines, and 70% in dash-dotted lines). A 1:1 mass ratio has been plotted for each panel in grey dotted lines.
In the text
thumbnail Fig. 9.

Surface rotation (top left), HRD (top right), surface 14N evolution (bottom left), and mass evolution (bottom right) of 50 M⊙ models at 10% Z⊙ for a range of critical rotation rates (0% in solid red lines, 40% in dash-dotted green lines, and 70% in dotted grey lines).
In the text
thumbnail Fig. 10.

Surface rotation (top left), HRD (top right), surface 14N evolution (bottom left), and mass evolution (bottom right) of 300 M⊙ models at 10% Z⊙ for a range of critical rotation rates (0% in solid red lines, 40% in dash-dotted green lines, and 70% in dotted grey lines).
In the text
thumbnail Fig. 11.

Surface evolution of [Na/Fe] and [O/Fe] for 80–500 M⊙ models during core H burning at 1/30th Z⊙. The red circles represent detections of stars in NGC 5904, while blue triangles represent upper limits in [O/Fe] as defined by Carretta et al. (2009). The green dashed circle establishes the abundance ratios of field stars ([O/Na] = 0.6), with stars outside of it representing P2 stars.
In the text
thumbnail Fig. 12.

Ratio of [Na/Fe] as a function of [O/Fe] for 80–500 M⊙ models during core He burning at 1/30th Z⊙. First and second population stars from the cluster NGC5904, adopted from Carretta et al. (2009), are shown in red circles and blue triangles, respectively.
In the text
thumbnail Fig. 13.

Time evolution of surface abundances of a 80 M⊙ model at 3% Z⊙ calculated from the onset of core H burning until core He exhaustion. The [Na/Fe] and [O/Fe] are included.
In the text
thumbnail Fig. 14.

Time evolution of surface abundances of a 400 M⊙ model at 3% Z⊙ calculated from the onset of core H burning until core He exhaustion. The [Na/Fe] and [O/Fe] are included.
In the text
thumbnail Fig. 15.

Cumulative ejected masses of [Na/Fe] and [O/Fe] for 80–500 M⊙ models calculated at 1/30th Z⊙. Observations of stars in NGC 5904 are shown as in Fig. 11, adopted from Carretta et al. (2009). The insert includes the wide range of O-enrichment ejected by the 400–500 M⊙ models during the post-MS.
In the text
thumbnail Fig. 16.

Diluted ejected from 80–500 M⊙ models calculated at 1/30th Z⊙. Observations of stars in NGC 5904 are shown as in Fig. 11, adopted from Carretta et al. (2009). Dilution factors ranging from 1 to 100 are shown for each model (assorted symbols) and connected for each f value (coloured lines).
In the text
thumbnail Fig. 17.

Diluted ejected from initially 40% critically rotating 80–300 M⊙ models calculated at 1/30th Z⊙. Observations of stars in NGC 5904 are shown as in Fig. 11, adopted from Carretta et al. (2009). Dilution factors ranging from 1 to 100 are shown for each model (assorted symbols) and connected for each f value (coloured lines).
In the text
thumbnail Fig. 18.

Cumulative ejected masses of [Mg/Fe] and [Al/Fe] for 80–500 M⊙ models calculated at 1/30th Z⊙. Observations of stars in M2 are shown by blue triangles. The [Mg/Fe] ratio includes isotopes of 24, 25, 26Mg and 26Al, while the [Al/Fe] ratio is calculated for 27Al.
In the text
thumbnail Fig. A.1.

Central temperature evolution for 100 M⊙ models during core H burning calculated at various Z.
In the text
thumbnail Fig. A.2.

Surface rotation (top left), HRD (top right), surface 14N evolution (bottom left), and mass evolution (bottom right) of 50 M⊙ models at ZSMC for a range of critical rotation rates (0% in solid red lines, 40% in dash-dotted green lines, and 70% in dotted grey lines).
In the text
thumbnail Fig. A.3.

Surface rotation (top left), HRD (top right), surface 14N evolution (bottom left), and mass evolution (bottom right) of 300 M⊙ models at ZSMC for a range of critical rotation rates (0% in solid red lines, 40% in dash-dotted green lines, and 70% in dotted grey lines).
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.