The effects of surface fossil magnetic fields on massive star evolution: V. Models at low metallicity Open Access

ABSTRACT

At metallicities lower than that of the Small Magellanic Cloud, it remains essentially unexplored how fossil magnetic fields, forming large-scale magnetospheres, could affect the evolution of massive stars, thereby impacting the fundamental building blocks of the early Universe. We extend our stellar evolution model grid with representative calculations of main-sequence, single-star models with initial masses of 20 and 60 M|$_\odot$|⁠, including appropriate changes for low-metallicity environments (⁠|$Z = 10^{-3}$|–|$10^{-6}$|⁠). We scrutinize the magnetic, rotational, and chemical properties of the models. When lowering the metallicity, the rotational velocities can become higher and the tendency towards quasi-chemically homogeneous evolution increases. While magnetic fields aim to prevent the development of this evolutionary channel, the weakening stellar winds lead to less efficient magnetic braking in our models. Since the stellar radius is almost constant during a blueward evolution caused by efficient chemical mixing, the surface magnetic field strength remains unchanged in some models. We find core masses at the terminal-age main sequence between 22 and 52 M|$_\odot$| for initially 60 M|$_\odot$| models. This large difference is due to the vastly different chemical and rotational evolution. We conclude that in order to explain chemical species and, in particular, high nitrogen abundances in the early Universe, the adopted stellar models need to be under scrutiny. The assumptions regarding wind physics, chemical mixing, and magnetic fields will strongly impact the model predictions.

1 INTRODUCTION

Understanding the physics of massive stars (⁠|$M_{\rm ini} \gt 8$| M|$_\odot$|⁠) is of paramount importance for studies of the early Universe (e.g. Garcia et al. 2021; Eldridge & Stanway 2022; Klessen & Glover 2023; Vink et al. 2023). The role of massive stars in cosmic re-ionization, chemical enrichment, and stellar feedback makes them crucial building blocks over various spatial and temporal scales (e.g. Chiaki, Susa & Hirano 2018; Chiaki & Wise 2019; Chiaki et al. 2020). In low-metallicity environments, from here on defined as metallicity lower than that of the Small Magellanic Cloud (SMC), that is, |$Z \lt 2.6\cdot 10^{-3}$| or |$Z/Z_\odot \lt 20~{{\ \rm per\ cent}}$| (e.g. Mokiem et al. 2007; Dopita et al. 2019), the detailed observational characterization of massive stars is still in its early stages. The first large sub-SMC catalogue of OB stars was recently compiled by Lorenzo et al. (2022) and a census of OBe stars in three low-metallicity dwarf galaxies was studied by Schootemeijer et al. (2022). Several other studies scrutinize individual or small number of metal-poor massive stars, nevertheless with important implications for physical processes (e.g. Ji et al. 2024; Mardini et al. 2024).

Currently, the source of high nitrogen content in distant galaxies – e.g. a factor of 3 nitrogen enrichment compared to the solar value in GN-z11 observed by the James Webb Space Telescope – remains elusive, with Wolf–Rayet type stars, very massive stars and supermassive stars suggested as their possible origin (e.g. Charbonnel et al. 2023; Vink 2023; Marques-Chaves et al. 2024; Watanabe et al. 2024). Nitrogen enrichment from stars is also reported for several other low-metallicity dwarf galaxies (e.g. Yarovova et al.2023).

Stellar evolution models are commonly utilized to interpret spectroscopic observations and guide galactic and cosmological-scale studies. Various codes and assumptions have been employed to calculate grids of low and zero-metallicity massive star evolutionary models (e.g. Meynet & Maeder 2002b; Heger & Woosley 2010; Takahashi, Umeda & Yoshida 2014; Szécsi et al. 2015; Limongi & Chieffi 2018; Groh et al. 2019; Martins & Palacios 2021; Murphy et al. 2021; Tsai et al. 2023; Roberti, Limongi & Chieffi 2024; Volpato et al. 2024). Insofar, some physical ingredients in massive star evolutionary models, namely, stellar winds and rotation, are often extrapolated from constraints available in Galactic environments. Moreover, until recently evolutionary models of massive stars – both at high and low metallicity – seldom considered surface magnetic fields. Therefore, the impact of this physical mechanism remains mostly without quantification. For Population III stars with zero metallicity, Haemmerlé & Meynet (2019) calculated supermassive star models with magnetic braking and Lou & Ma (2022) investigated magnetic equilibria and radial pulsations. Urpin (2019) studied the thermal induction of magnetic fields. For low metallicities of |$Z/Z_\odot \sim 1/10$| and |$1/50$|⁠, Li et al. (2023) utilized stellar evolution calculations, including a branch of magnetic models with a 1 kG field strength, to study supernova progenitors. Nevertheless, until now, a coherent extension of typical massive star models with surface magnetic field effects at sub-SMC metallicities is missing. This underpins a major uncertainty in evolutionary model calculations since both stellar winds and rotation are drastically impacted by large-scale magnetic fields, which form a magnetosphere around the star (e.g. ud-Doula & Owocki 2002; Townsend & Owocki 2005; ud-Doula, Owocki & Townsend 2009; Owocki et al. 2016).

This paper is part of a series in which we aim to explore the effects of surface fossil magnetic fields on massive star evolution (Keszthelyi et al. 2019, 2020, 2021, 2022, Papers I–IV from here on). Papers I, II, and III focused on solar metallicity models, while in Paper IV we computed a large grid for metallicities representative of the solar neighbourhood and the Magellanic Clouds. For a very brief summary of the series, we refer the reader to the introduction of Paper IV. Given that stellar magnetism is known to largely impact the evolution of stars at high-metallicity environments, it is natural to wonder how stellar model predictions at low and zero metallicity are impacted by magnetism. In this work, we utilize our implementation of surface fossil magnetic field effects to study typical massive star evolutionary models with initial masses of 20 and 60 M|$_\odot$| at a domain of |$Z = 10^{-3}{\!-\!}10^{-6}$|⁠, corresponding to |$Z/Z_\odot \approx 6\,\times \, 10^{-2}{\!-\!}6 \,\times \, 10^{-5}$|⁠. In this metallicity range, the incidence rate and general magnetic properties of massive stars are not yet known observationally since spectropolarimetric observations remain limited to detections in a Galactic environment (e.g. Bagnulo et al. 2020). In addition to the parametrization and schemes introduced in Paper IV, we consider appropriate changes for low-metallicity models. This includes two branches to calculate the applied mass-loss rates. This is because stellar wind physics becomes more uncertain and observationally unconstrained in sub-SMC metallicities. In fact, when lowering the metallicity a critical limit may be reached, leading to effectively no mass-loss. We discuss these modelling ingredients in Section 2.

The paper is organized as follows. In Section 2, we detail the main assumptions and general model set-up. In Sections 3, 4, and 5, we present and scrutinize the stellar structure and evolution models from our computations. In Section 6, we discuss the implications, and in Section 7, we outline future avenues. Then, we conclude our findings in Section 8.

2 METHODS

2.1 General model set-up

We use the software instrument Modules for Experiments in Stellar Astrophysics mesa release 22.11.01 (Paxton et al. 2011, 2013, 2015, 2018, 2019) and Software Development Kit (sdk) version 22.6.1 (Townsend 2020). The mesa microphysics are summarized in the Appendix. The models computed in this work are single-star models covering the main-sequence evolution. Table 1 lists the main parameter space. Here, we summarize some main points of the model set-up but refer the reader to Paper IV for further details.

2.1.1 Initial mass

In this study, we consider initial masses of |$M_{\rm ini} = \lbrace 20, 60 \rbrace \, \mathrm{M}_{\odot }$|⁠. We choose this mass range for four reasons. (i) It is consistent with the upper mass range presented in Paper IV, where the effects of winds and rotation are the most pronounced. (ii) Stable fossil fields are anchored in the radiative envelope and are most likely expelled from the convective cores of main-sequence massive stars (Keszthelyi 2023). For initial masses above 60 M|$_\odot$|⁠, the convective core starts to occupy a much larger fraction of the total mass of the star. Due to strong stellar winds that remove much of the radiative envelope, the core boundary will reside much closer to the stellar photosphere. In addition, the envelopes of Very Massive Stars (VMSs) may reach super-Eddington luminosities. In those layers, convective turbulence could develop to transport the luminosity (Quataert & Shiode 2012; Quataert et al. 2016,¹). The stability of fossil fields is unclear under those conditions (see Section 7). (iii) Given the initial effective temperatures (⁠|$\approx 50$| kK), the models are most representative of typical O-type stars. Therefore, these models may be directly compared with observations of stars in metal-poor galaxies. (iv) VMSs and supermassive stars (SMSs) might dominate the luminosity budget and chemical evolution of stellar populations and have received considerable interest at near zero-metallicity conditions. The ‘typical’ massive stars that we model in this work would be just as frequent in number as more exotic ones if the initial stellar mass function (IMF) was completely flat for high masses. However, the IMF is likely not flat but top-heavy close to zero metallicity, supported by star formation studies (e.g. Chiaki & Yoshida 2022; Chon et al. 2022, see also Bastian, Covey & Meyer 2010). In addition, protostellar feedback could also play a role in setting the upper mass limit of a forming star (Hosokawa et al. 2011). These imply that typical massive stars are possibly more numerous than VMSs and SMSs and therefore could have an important contribution to the luminosity budget and chemical evolution.

2.1.2 Initial magnetic field strength

We consider initial equatorial surface magnetic field strengths |$B_{\rm eq,~ini} = \lbrace 0, 0.5, 1, 10, 50 \rbrace $| kG. The upper limit is the same as in Paper IV and is chosen as an exploration of the parameter space. While observationally inferred magnetic field strengths of Galactic massive stars are in the kG range (Grunhut et al. 2017; Shultz et al. 2018), it is presently unclear if fossil fields in low-metallicity environments would have different strengths. Star-forming molecular clouds have more than sufficient magnetic flux to produce magnetized stars (e.g. Keszthelyi 2023, and references therein). At high metallicity, efficient dissipation via non-ideal magnetohydrodynamic processes need to be invoked to reduce the magnetic flux from order of |$10^{33}$| G cm|$^{2}$| in star-forming cores to around |$10^{27}$| G cm|$^{2}$| that is observed in stars (Petit et al. 2017, Paper II). So far, it remains unclear what initial magnetic flux and dissipation could characterize molecular clouds at low metallicities.

2.1.3 Initial rotation

Consistently with Paper IV, rotation is initialized assuming a flat angular velocity profile and we use the mesa control option |$\Omega /\Omega _{\rm crit}$| (set to 0.5), where the critical angular velocity is defined in mesa as |$\Omega _{\rm crit}~=~\sqrt{(1-\Gamma) GM/R }$| with |$\Gamma$| the Eddington parameter calculated for all opacity sources at the stellar photosphere. This results in initial surface equatorial rotational velocities in the range of |$\approx$| 350–550 km s|$^{-1}$|⁠, with increasing values towards higher initial masses and lower metallicity (due to smaller stellar radii). This assumption is a reasonable approximation since star formation studies indicate that low-metallicity massive stars may initiate their evolution with high rotational velocities (Hirano & Bromm 2018). We refrain from computing models with various initial |$\Omega /\Omega _{\rm crit}$| ratios since to some extent this is substituted by our assumptions regarding chemical mixing.

2.1.4 Chemical mixing schemes

As in Paper IV, we adopt two chemical mixing schemes that account for transporting elements inside the stellar model. The ‘Mix1’ scheme contains the typical description used in mesa, which is based on the works of Endal & Sofia (1978) and Pinsonneault et al. (1989). The ‘Mix2’ scheme is our implementation of the Zahn (1992) formalism into mesa. We found the latter case to be very efficient in chemical mixing at higher metallicity environments, therefore these scenarios may be considered as boundary cases for weak and strong chemical mixing. We further elaborate on this, particularly, in Section 3.

2.1.5 Angular momentum transport and braking schemes

In Paper IV, the models were split into 3 main branches to account for uncertainties in applying angular momentum transport and magnetic braking. ‘INT’ models are assumed to be solid-body rotating. Angular momentum loss accounting for magnetic braking and stellar winds is distributed through all layers of these models. ‘SURF’ models on the other hand have a rotation profile, which allows for radial differential rotation. In the SURF models, specific angular momentum is removed from the envelope and surface layers in one time-step. It is assumed that this region contains 20 per cent of the total mass of the star. In the ‘NOMAG’ models, we assume that only hydrodynamic instabilities transport angular momentum in mesa’s diffusive scheme. In this case, angular momentum loss follows the default mesa procedure. Specific angular momentum is subtracted from those near surface layers where mass is removed to account for stellar winds.

2.1.6 Differences compared to Paper IV

While the models are generally consistent with those of Paper IV, we should highlight that in the present work we use a newer mesa release and consider a more extended 21-isotope ‘approx21.net’ nuclear network compared to the 8-isotope ‘basic.net’ used in Paper IV.

2.2 Main assumptions for low metallicity

To consider models at low metallicity, we rely on some further assumptions. We now discuss changing the initial chemical composition, as well as the consideration of magnetic braking and stellar-wind-driven mass-loss.

2.2.1 Metallicity

We adopt initial metallicities in our models from |$Z = 10^{-3}$| to |$10^{-6}$|⁠, which is approximately from |$\approx$|1/17 to |$\approx$|/17 000 of solar metallicity, considering |$Z_{\odot } \approx 0.014{\!-\!}0.017$| (Asplund et al. 2009; Magg et al. 2022). This range encompasses the presently known most metal-poor massive stars observed in Sextans A and IZw18 and also extends to the high-redshift universe. To obtain the corresponding initial hydrogen (⁠|$X_{\rm ini}$|⁠) and initial helium (⁠|$Y_{\rm ini}$|⁠) abundances, we test a few different interpolation methods. To this extent, we use anchoring values of metals and helium of the primordial abundance, the initial abundance in the Large and Small Magellanic Clouds, and the initial solar abundance. Table 2 lists the mean values, where the study of Dopita et al. (2019) used the mean values of nine different investigations and the study of Magg et al. (2022) considered the mean of six values from different studies and methods. Generally, the helium abundances at a given metallicity are consistent between various studies within a 5 per cent level.

A usual method is to perform a two-point linear interpolation, typically using the proto-solar and primordial helium abundances. We additionally consider whether four-point linear and cubic interpolations (including data from the Magellanic Clouds) could result in differences and find that in the metallicity range of our interest, the three methods agree very well with deviations below a 0.1 per cent level (Fig. 1). We opt to use the results obtained from the four-point linear interpolation.

The resulting initial helium and hydrogen (⁠|$X = 1-Z - Y$|⁠) abundances are listed in Table 3. We assume that all initial hydrogen is in the |$^{1}$|H isotope. For simplicity, we make the assumption that the initial helium abundance is distributed between its two isotopes as |$^{3}\mathrm{He}/^{4}\mathrm{He} \sim 10^{-4}$| across all metallicities (see e.g. Meynet & Maeder 2002b). We tested models with initially only |$^{4}\mathrm{He}$| and found that they generally lead to small differences in the main model characteristics. None the less, such an assumption can indeed be important when extrapolating to zero metallicity environments (c.f. West, Heger & Austin 2013; Takahashi 2018). Although it is likely that the distribution of individual metal abundances over cosmic times does not follow a simple scaling down from measured solar values, the uncertainty in describing individual metal elemental abundances makes it impossible to tailor this input in our models. For this reason, we choose to use the initial distribution of metals with the option initial_zfracs = 8 in mesa, corresponding to the Asplund et al. (2009) data, modified by Nieva & Przybilla (2012) and Przybilla et al. (2013), considering B-type stars in the solar neighbourhood.² In all cases, we use the Lodders (2003) isotopic ratios and the metallicity adopted for the chemical evolution is the same value that is used for the opacity calculations.

2.2.2 Magnetic braking

Large-scale magnetic fields, which form a magnetosphere around the star, are efficient in transporting angular momentum along the field lines independently of the plasma flow. As the star rotates, this Maxwell stress and associated Poynting flux can transfer angular momentum from the magnetic field to the surrounding wind material. The angular momentum loss from the star leads to rapid spin down. This is well evidenced by the generally slow surface rotation of magnetic stars (e.g. Shultz et al. 2018, 2019).

Two-dimensional magnetohydrodynamic simulations of rotation-aligned dipole fields of massive stars by ud-Doula et al. (2009) indicated that the angular momentum carried out by a magnetically torqued stellar wind follows the same simple, split-monopole scaling law derived for the Sun by Weber & Davis (1967), |$\dot{J} =\frac{2}{3} {\dot{M}} \, \Omega \, R_A^{2}$| – with, however, the Alfvén radius |$R_A$| now given by the dipole scaling |$R_A \sim \eta _\star ^{1/4}$|⁠, instead the oft-quoted, stronger scaling |$R_{A} \sim \eta _\star ^{1/2}$| for a split monopole. Here, |$\eta _\star$| is the ‘wind magnetic confinement parameter’ and represents the ratio of magnetic energy density and wind kinetic energy density (ud-Doula & Owocki 2002). For typical confinement by a dipole field (⁠|$\eta _\star \approx 10 - 100$|⁠), |$\eta _\star ^{1/2} \propto B \dot{M}^{-1/2} \propto R_A^2$|⁠, with B the magnetic field strength and |$\dot{M}$| the mass-loss rate. Therefore, the angular momentum loss |$\dot{J}$| is proportional to |$B \dot{M}^{1/2}$|⁠. This scaling forms the basis for modelling the rotational evolution of magnetic massive stars (Paper I, Paper II, Paper III, Paper IV).

Thus far, it is little explored how magnetic braking operates in low wind density conditions. The key point to achieve the stellar spin-down is that the magnetic field at the distance of the Alfvén radius has to be able to transfer angular momentum to some surrounding material that is detached from the star. In principle, this could also be the interstellar medium. In our models, we assume that magnetic braking operates for (very) low mass-loss rates; however, we do not apply it in the case of |$\dot{M}=0$|⁠.

2.2.3 Stellar wind physics

Radiative-line driving is a well-understood mechanism that leads to significant mass-loss from high-metallicity, hot, massive stars (e.g. Puls, Vink & Najarro 2008; Vink 2022, for reviews). However, the nature of stellar winds at metallicities lower than that of the SMC is highly uncertain. Radiative-line driving is expected to become less efficient when lowering the metallicity, primarily due to the decrease in the iron abundance, which is the major driver at near-solar metallicity conditions.

At low wind densities (as predicted for low metallicity environments and/or for stars of comparatively low luminosity), a potential de-coupling of ions from the bulk wind plasma could lead to a runaway effect. The corresponding frictional heating might dominate the energy balance, thus even enforcing the de-coupling (Springmann & Pauldrach 1992). In a series of publication, Krtička & Kubát (2000, 2001a, b, 2002) and Krtička, Kubát & Groote (2006) studied two- and three-component stellar wind models (partly including frictional heating), and surprisingly found steady-state solutions even for low wind densities. They argued that the de-coupling of the accelerated ions from the bulk H+He plasma becomes quite difficult, if not impossible, under typical conditions. We note here that in contrast to metals, hydrogen, and helium have only a few lines available; therefore, any form of wind driving from H+He alone would result in very low mass-loss rates (Krtička & Kubát 2006).

Owocki & Puls (2002) confirmed the findings by Krtička & Kubát (2002) for stationary assumptions, and discovered a new, very strong and rapidly propagating ‘ion separation instability’. This could destroy any stationary flow at those densities (or even somewhat higher ones) where the de-coupling may occur. Based on these results, we follow the parametrization of Springmann & Pauldrach (1992) to calculate the critical mass-loss rate |$\dot{M}_{\rm crit}$| as a function of metallicity where the metal ions are predicted to de-couple from the H/He bulk wind. Below this limit, the accelerated metal ions might de-couple from the bulk wind, such that mass-loss either terminates at all, or continues at extremely low rates (if driving by H and He alone was possible). We have also compared with the alternative approach by Owocki & Puls (2002), which results in slightly different values, though at maximum by a factor of two. For simplicity, and to obtain a lower limit for |$\dot{M}_{\rm crit}$|⁠, we neglect the impact of frictional heating.

We start with the condition for the boundary between coupling and de-coupling, by equating

where |$\Gamma _{\rm B}$| is the |$\Gamma$|-value (⁠|$= g_{\rm rad}(r)/g_{\rm grav}(r)$|⁠) for which the drift between ions and bulk plasma reaches the thermal speed. For larger drifts, the collisional coupling decreases rapidly, and de-coupling becomes (principally) possible in a runaway fashion. |$\Gamma _{\rm L}$| is the radiative line force scaled to the total wind density, divided by the gravitational acceleration. In the supersonic approximation and neglecting the acceleration by electrons (this is one of the differences compared to the study by Owocki & Puls 2002),

Following equations (12) and (13) of Springmann & Pauldrach (1992), i.e. neglecting the number density of metals compared to the number density of the bulk plasma (with a Helium content).³|$Y_{\rm He}=n_{\rm He}/n_{\rm H}$|⁠, and a number of |$I_{\rm He}$| free electrons per He nucleus), |$\Gamma _{\rm B}$| can be approximated by

with |$Y_i$| the metallic abundance (expressed as a number fraction, |$Y_i=n_i/n_{\rm H}$|⁠), |$\ln \Lambda$| the Coulomb logarithm (⁠|$\approx 20)$|⁠, |$Z_i$| the average charge of the metallic ions in units of the electron charge e (⁠|$Z_i \approx$| between 2 and 3 for OB stars), |$m_{\rm H}$| the hydrogen mass, |$k_{\rm B}$| the Boltzmann constant, |$G_{\rm max} = 0.214$| the maximum value of the Chandrasekhar function, and all other quantities as conventionally designated.

The velocities and velocity gradients can be expressed by the canonical |$\beta$| law (e.g. Pauldrach, Puls & Kudritzki 1986), which for the prototypical value |$\beta = 1$| results in

with |$v_{\infty }$| the terminal speed of the (coupled) wind, and |$R_\star$| the stellar radius. Equating the relations for |$\Gamma _{\rm B}$| and |$\Gamma _{\rm L}$|⁠, and requiring a completely de-coupled wind (i.e. setting |$v(r) =: v_{\rm crit} \approx 0.1 v_{\infty }$| as a prototypical value for the speed of a line-driven wind at its critical point, see Castor, Abbott & Klein 1975; Pauldrach et al. 1986), we finally obtain

Overall, the critical mass-loss rate (w.r.t. de-coupling) scales as

where the temperature might be replaced by the effective one. Obviously, the lower the metallicity, the larger |$\dot{M}_{\rm crit}$|⁠. Accounting for the fact that (observed) terminal velocities of OB stars scale with |$T_{\rm eff}$| (Hawcroft et al. 2023), the (effective) temperature has a large impact on the critical mass-loss rate.

Our results for a typical 60 M|$_\odot$| star near the ZAMS are shown in Fig. 2. If the actual mass-loss rate is below the de-coupling limit, then the radiative force might be insufficient to drive a wind. Here, the actual mass-loss rate is assumed with a typical value for Galactic massive stars and then scaled down by a 0.85 exponent for lower metallicity (Vink, de Koter & Lamers 2001; Mokiem et al. 2007). This results in our first-order estimate shown with the magenta line.

From solar metallicity to a metallicity of |$Z/Z_\odot \sim 6 \,\times \, 10^{-3}$|⁠, the critical mass-loss rates for O-type dwarfs with effective temperatures of 60 kK are on the order of |$10^{-11}$| and |$10^{-8}~$| M|$_\odot$| yr|$^{-1}$|⁠, respectively. Typical mass-loss rates calculated from wind prescriptions or obtained from observational diagnostics are above this limit at solar metallicity. However, at around |$Z/Z_\odot \sim 6 \,\times \, 10^{-3}$| and at even lower metallicities, the calculated rates will tend to fall below the de-coupling limit. Thus, according to our estimates, reaching the de-coupling limit will take place in the metallicity range considered in our study. For this reason, we implement the calculation of the critical de-coupling limit in the mesarun_star_extras extension and systematically compare whether the mass-loss rate calculated from a prescription is above or below this limit at each point during the evolution of a stellar model. In the following, we use the terminologies of calculated and applied mass-loss rates. By calculated rates, we refer to theoretical predictions from a wind prescription. By applied rates, we refer to the rates used to reduce the mass of the stellar model.

2.2.4 Mass-loss schemes

Given the above considerations, we distinguish between two modelling scenarios. In the first case (W1), even if the calculated mass-loss rates are below the critical value, we assume that some other form of wind-driving might take place, thus we apply non-zero mass-loss. Even if radiative-line driving becomes completely inefficient, winds from low and zero metallicity massive stars might be launched by other mechanisms (which mechanisms might not depend on metallicity). The prime candidates are pulsations-driven mass-loss (e.g. Appenzeller 1970; Shiode, Quataert & Arras 2012), Super-Eddington winds (e.g. van Marle, Owocki & Shaviv 2008; Sanyal et al. 2015; Owocki, Townsend & Quataert 2017), or other non-steady state eruptive mass-loss episodes. Though, these mechanisms might require specific conditions (evolutionary stage, high luminosity, etc.). For the sake of simplicity, we apply the originally calculated mass-loss rates in the W1 branch. In the second case (W2), we assume that the models experience no mass-loss at all if their calculated rates would be below the de-coupling limit by applying |$\dot{M}= 0$|⁠. Note however that this is time dependent and may not characterize the entire evolution of a model.

Rapid rotation could enhance (the existing) mass-loss. We include this effect and calculate rotational enhancement on the mass-loss rates following Maeder & Meynet (2000). In the case of models approaching critical rotation (which threshold we arbitrarily set to |$\Omega /\Omega _{\rm crit} = 0.9$|⁠), we apply a numerical scheme built into mesa to artificially increase the mass-loss rates. Insofar, the application of such numerical techniques is seldom verified and we will refrain from discussing physical details near critical rotation, which is a largely unexplored and unresolved issue.

In consistency with Paper IV, we adopt the Vink et al. (2001) rates reduced by a factor of 2 for hot, hydrogen-rich stars (⁠|$T_{\rm eff} \gt 10$| kK, |$X\gt 0.4$|⁠), and the Nugis & Lamers (2002) rates for hot, helium-rich stars (⁠|$T_{\rm eff} \gt 10$| kK, |$X\lt 0.4$|⁠). If chemical mixing is very efficient in the models (or the initial rotation is high), they can become helium rich and meet this condition already on the main sequence. We tested the new Sander & Vink (2020) rates (developed for |$2.0 \gt Z/Z_{\odot } \gt 0.02$|⁠); however, the exponential decrease in the predicted mass-loss rates below SMC metallicity leads to prohibiting numerical issues. We use the metallicity scaling included in the Vink et al. (2001) prescription, that is |$\dot{M} \sim (Z/Z_\odot)^{0.85}$|⁠. In this case, we use the initial metallicity of the models. The Nugis & Lamers (2002) prescription also includes a metallicity scaling (⁠|$\dot{M} \sim Z(t)^{0.5}$|⁠), which however uses the actual (time-dependent) metallicity of the stellar model (see discussion by Higgins et al. 2021).

3 DEGENERACY BETWEEN CHEMICAL MIXING EFFICIENCY AND INITIAL ROTATIONAL VELOCITY

In this section, we briefly demonstrate how varying either chemical mixing or the rotational velocity can lead to similar effects. Traditionally, uncertainties in chemical mixing have been less explored in evolutionary models (however, see e.g. Meynet et al. 2013; Nandal et al. 2024), instead the initial rotation rates were used to study various scenarios. These findings indicated that when rotation is rapid enough, the star may have a structure that is almost perfectly mixed in composition (Maeder 1987; Maeder & Zahn 1998; Yoon, Langer & Norman 2006; Cantiello et al. 2007). Such quasi-chemically homogeneous evolution (QCHE) is expected to be more common at metallicities lower than that of the SMC (e.g. Langer 2012; Kubátová et al. 2019). In Paper IV, we fixed the initial rotation rates and tested two cases for chemical mixing, the less efficient Mix1 scheme and the more efficient Mix2 scheme.

In Fig. 3, we demonstrate that the Mix1 case tends toward the Mix2 case if the initial angular velocity is increased. The models considered here are without magnetic field effects and in the W1 wind scheme (mass-loss rates are applied even if the de-coupling limit is reached). Keeping the initial rotation the same (⁠|$\Omega /\Omega _{\rm crit}=0.5$|⁠), the Mix1 model evolves redward, whereas the Mix2 model evolves blueward in the HRD (top left panel). The latter evolutionary trajectory is indeed associated with QCHE. A model within the Mix1 chemical mixing scheme but with higher initial rotation (⁠|$\Omega /\Omega _{\rm crit}=0.8$|⁠) has a Zero Age Main Sequence (ZAMS) that is slightly shifted to a lower effective temperature. This is because the centrifugal force adds support against gravity and allows for a larger stellar radius. This track (orange dashed line) follows an evolution that is almost identical to the Mix2 model which is initially less rapidly rotating but has a higher efficiency of transporting chemical species. Besides the obvious quantitative difference in rotational velocities, the evolution of the surface abundances of He and CNO elements in the Mix1 |$\Omega /\Omega _{\rm crit}=0.8$| model also tends toward the Mix2 |$\Omega /\Omega _{\rm crit}=0.5$| case. For this reason, in the following it may be considered that the Mix2 scheme with initially |$\Omega /\Omega _{\rm crit}=0.5$| gives an indication for results where, instead of the mixing efficiency, the rotational velocity was increased in a Mix1 model. With our approach, we wish to emphasize that chemical mixing is largely unresolved and as such remains a key physical aspect to calibrate in evolutionary models.

4 RESULTS WITHIN THE W1 WIND SCHEME

In this section, we present our findings when utilizing the assumption of the W1 wind scheme. This results in gradually weaker stellar winds as a function of decreasing metallicity. In the next section, we investigate the impact of this choice by discussing the W2 wind scheme, where the applied mass-loss rates can become zero.

4.1 Evolution of main stellar parameters and surface abundances

In contrast to high-metallicity environments, stellar mass-loss is modest, whereas rotation and chemical mixing lead to much more pronounced effects at low metallicity (e.g. Meynet & Maeder 2002a). We proceed with discussing our models with the assumption of a fixed initial rotation rate and two boundary cases (Mix1, Mix2) for chemical mixing.

Fig. 4 shows the HRD for models with an initial mass of 60 M|$_\odot$| at the initial metallicities considered in our study (⁠|$Z=10^{-3}$|–|$10^{-6}$|⁠, from top left to lower right panels), using the W1 wind scheme. The initially 20 M|$_\odot$| models follow similar qualitative trends. Here, we only focus on the magnetic models in the INT scheme (solid-body rotators). The SURF models shows similar evolutionary trajectories. Following Szécsi et al. (2015), we may distinguish between three evolutionary channels: (i) a ‘classical’, redward trajectory in the HRD, (ii) a transitional, intermediate situation, in which some part of the evolution is blueward but the tracks return to a redward trajectory, and (iii) a completely blueward evolution. The blueward evolution may be associated with QCHE, that is, the star is almost completely mixed in composition. In the following, we simply use the initial, ZAMS effective temperature of the models as a reference and compare whether the tracks follow an evolution towards lower (redward) or higher (blueward) effective temperatures on the main sequence.

The initial metallicity plays a large role in the development of the QCHE channel. Lowering the metallicity leads to the well-known effect of shifting the ZAMS effective temperatures to higher values, which is due to the smaller stellar radius. For this reason, and consequently due to the higher rotational velocities, lower metallicity also favours a more extended blueward evolution towards high effective temperatures and luminosities.

Within the Mix1 chemical mixing scheme, the models tend to evolve along a classical path. However, when lowering the metallicity, the evolutionary trajectory can remain vertical in the HRD for up to half of the main sequence lifetime (⁠|$\sim 2$| Myr), or even with a slightly increasing |$T_{\rm eff}$|⁠. In the Mix1 case, the assumption about initial magnetic field strength or angular momentum transport scheme lead to relatively small differences in the evolutionary tracks, although the exact width of the main sequence, the lifetimes, and other physical quantities vary in the models (c.f. Papers II and IV).

Models with the more efficient Mix2 chemical mixing scheme tend to evolve bluewards in the HRD. The stronger the magnetic field, the slower the rotation becomes, resulting in less efficient mixing. Therefore, magnetic fields work to prevent the QCHE channel (c.f. Paper IV). Compared to a non-magnetic track, the magnetic tracks start to deviate from the QCHE channel and turn towards a classical evolution. However, while in higher metallicity environments weaker magnetic fields are sufficient to brake the rotation and suppress chemical mixing to prevent this evolutionary channel, in low-metallicity, a stronger magnetic field is required. In particular, at |${\it Z}=10^{-3}$| (top left) only the 10 and 50 kG models return to lower effective temperatures than the ZAMS value and follow a classical evolution. The model with 1 kG field strength follows a transitional evolution, initially evolving bluewards and returning to a redward trajectory towards the TAMS. At |$Z=10^{-6}$|⁠, only the 50 kG model in the INT angular momentum transport and magnetic braking scheme departs from the QCHE channel. This highlights that at lower metallicity environments even strong magnetic fields are insufficient to prevent the efficient chemical mixing.

The surface abundances of He and CNO elements are shown in Fig. 5 for the models with initial mass of 60 M|$_\odot$| at |${\it Z}=10^{-5}$|⁠. QCHE leads to helium-rich stars on the main sequence in the Mix2 case, except the 50 kG model which brakes its rotation and returns to a classical redward evolution (top left panel, c.f. also Fig. 4 and discussion above). In both Mix1 and Mix2 models, the stronger the magnetic field, the less helium enrichment is predicted. During the main sequence, as set by the CNO cycle, first the surface carbon and then the oxygen abundances decrease, while the surface nitrogen abundance increases rapidly. These changes are more prominent in the initially 60 M|$_\odot$| models, but the initially 20 M|$_\odot$| models also show similar trends. The magnetic models systematically deviate from the non-magnetic models. Stronger magnetic fields lead to less surface enrichment of nitrogen, whereas the depletion of surface carbon and oxygen is more moderate. This is because magnetic braking leads to slower spin, which decreases the efficiency of chemical mixing inside the star. We find that typical, order of kG magnetic fields tend to moderately impact the evolution of surface abundances, within about 0.2 dex. Very strong magnetic fields (50 kG) are needed to significantly inhibit chemical mixing.

4.2 Magnetic field evolution

Utilizing the assumption of magnetic flux conservation, the evolution of surface equatorial magnetic field strength only depends on the stellar radius (see Papers II and III and references therein). In Fig. 6, we compare the time evolution of the surface equatorial magnetic field strength of models with an initial mass of 60 M|$_\odot$| at the metallicities considered in our study. At |$Z = 10^{-3}$|⁠, all models within this specific setup follow a decline in their surface equatorial magnetic field strength. At |$Z = 10^{-4}$|⁠, the magnetic field evolution of the models with initial field strength of 0.5 and 1 kG in the Mix2-SURF and Mix2-INT schemes remains practically constant. At |${\it Z} = 10^{-5}$| and |$Z = 10^{-6}$|⁠, also some models with initially 10 and 50 kG field strength display a nearly constant magnetic field evolution in the Mix2 chemical mixing scheme.

This might seem surprising; however, the reason for these trends is the competition between chemical mixing and surface magnetism. Strong chemical mixing favours blueward evolution in the HRD with nearly no change in the stellar radius, whereas a strong magnetic field favours a classical evolution with increasing radius due to slow rotation and less efficient mixing. Consequently, our results have interesting implications. (i) At low metallicity, if chemical mixing is efficient (or initial rotation is rapid, c.f. Section 3), magnetic fields could maintain the same strength over the course of their entire main sequence evolution. This is unlike in observed Galactic samples, where a decline in magnetic field strength over the main sequence has been identified (Landstreet et al. 2007, 2008; Fossati et al. 2015; Shultz et al. 2019). Approximately, an order of magnitude decline in the surface equatorial magnetic field strength is also expected from our low-metallicity model predictions in the Mix1 scheme with less efficient chemical mixing. (ii) High-mass, main-sequence stars are likely to become Wolf–Rayet type stars during the post-main sequence. In fact, several of our models are also Wolf–Rayet type-like already on the main sequence in terms of their surface helium enrichment as shown in Fig. 5. Therefore, the initial equatorial magnetic field strengths of 0.5, 1, and 10 kG (implying a 1, 2, and 20 kG dipole field strength), which are predicted to remain quasi-constant during the main-sequence evolution in the Mix2 models at |$Z=10^{-5}$| and |$Z=10^{-6}$| (lower row of Fig. 6), might be assumed to be the magnetic field strength of a Wolf–Rayet type star inherited from its progenitor.

4.3 Evolution of surface rotation

Stellar rotation plays a major role in the evolution of low-metallicity massive stars. In diffusive angular momentum schemes as in mesa, the speed of rotation is the main driver of chemical mixing (primarily via meridional currents) in radiative envelopes.⁴ Therefore, it is important to quantify how magnetic fields impact rotation.

For this reason, we scrutinize the evolution of the surface equatorial rotational velocity as predicted in our models (Fig. 7). If the star has a magnetosphere, most of the angular momentum loss is driven by the magnetic field (Weber & Davis 1967; ud-Doula et al. 2009). At |$Z = 10^{-3}$|⁠, the non-magnetic Mix1 models maintain the highest rotational velocity on the main sequence. The non-magnetic Mix2 model follows a similar rotational evolution until about 3.5 Myr. At that point, the increasing surface helium abundance leads to Wolf–Rayet type mass-loss rates. These rates are much higher than the optically-thin wind prescriptions, therefore an enhanced rotational braking proceeds. In general, the Mix1 models can inflate their envelopes more, which also contributes to lowering their surface rotation due to angular momentum conservation. The Mix2 models tend to remain more compact, which favours maintaining higher surface rotational velocities. In some cases, the models can experience a spin up. This is because the time-scale to transport and replenish angular momentum in the surface layers is shorter than the time-scale to remove it. The magnetic models follow the trend evidenced in high-metallicity, namely, the stronger the magnetic field, the faster the surface rotation brakes (e.g. Papers II and IV). At |$Z = 10^{-5}$| and |$Z = 10^{-6}$|⁠, the rotational evolution becomes more complex due to the tendency towards QCHE. In these cases too, the sudden drop (by hundreds of km s|$^{-1}$|⁠) in surface rotational velocity close to the end of the main sequence is attributed to switching to Wolf–Rayet type mass-loss rates. Overall, we can see that when lowering the metallicity, progressively stronger magnetic fields are needed to brake rotation. The difference between the INT and SURF magnetic braking schemes is the most prominent for strong magnetic fields. For 10 and 50 kG fields, the INT models display lower surface rotational velocities over their evolution than the SURF models. This is because efficient chemical mixing remains in the SURF models, which allows them to evolve on a blueward evolutionary trajectory. In turn, these models are more compact and therefore can maintain higher surface rotational velocities than the INT models.

5 RESULTS WITHIN THE W2 WIND SCHEME

In this section, we use a few selected examples to highlight how the assumption on stellar wind impacts the results. We begin with discussing the scenarios, in which the de-coupling limit is reached and zero mass-loss rates are applied in the W2 scheme.

5.1 Mass-loss rates

As shown in Section 2.2.3, the de-coupling limit concerns the metallicities of |$Z = 10^{-5}$| and |$Z = 10^{-6}$| in the case of an initial mass of 60 M|$_\odot$|⁠. When decreasing the initial mass, the lower luminosity results in weaker winds. Consequently, the calculated mass-loss rate of the 20 M|$_\odot$| model falls below the de-coupling limit already at |$Z = 10^{-4}$|⁠. However, the critical mass-loss rate to reach the de-coupling of ions and bulk plasma in the stellar wind is dependent on stellar parameters, which continuously change during the evolution of the star. In Fig. 8, we demonstrate this with the specific example of models using the Mix1-INT scheme with |$B_{\rm eq,~ini}= 1$| kG. At |$Z = 10^{-4}$|⁠, the model with |$M_{\rm ini}=20$| M|$_\odot$| initiates its evolution with zero applied mass-loss rates since the calculated rates (orange line) are below the critical de-coupling limit (grey line). At around 7 Myr, the calculated rates reach above the de-coupling limit hence the model will experience mass-loss from that time onwards. At |$Z = 10^{-5}$|⁠, the calculated mass-loss rate for the |$M_{\rm ini}=20$| M|$_\odot$| model is constantly below the de-coupling limit, thus we apply zero mass-loss in that case. Contrary to this, at |$Z = 10^{-4}$|⁠, the calculated mass-loss rate for the |$M_{\rm ini}=60$| M|$_\odot$| model is above the de-coupling limit, hence mass-loss is applied. In this case, the models in the W1 and W2 branches are completely identical. At |$Z = 10^{-5}$|⁠, the |$M_{\rm ini}=60$| M|$_\odot$| model in the W2 has zero applied mass-loss rates since the calculated values are systematically below the critical de-coupling limit throughout the main-sequence evolution.

5.2 Critical rotation

One of the immediate consequences of applying zero mass-loss rates is that our massive star evolutionary models predict a spin-up. Therefore, these models are much more prone to approach the limit of critical rotation, at which point the centrifugal force balances gravity and radiation pressure (Maeder & Meynet 2000).

The physics of critical rotation and its implementation in one-dimensional stellar evolution calculations have remained challenging, unresolved problems. Here, we simply aim to identify the parameter space that favours the development of critical rotation, but we refrain from further discussion of this complex issue. To this extent, Figs 9 and 10 show the tendency towards critical rotation for initially 60 M|$_\odot$| models using the W1 and W2 wind schemes. None of the models reach critical rotation when using the W1 scheme. The mass-loss rates are enhanced when reaching |$\Omega /\Omega _{\rm crit} \approx 0.9$| and this prevents reaching critical rotation. In general, lower metallicity and weaker magnetic fields favour the development of near critical rotation. At |$Z = 10^{-3}$| and |$Z = 10^{-4}$|⁠, the W1 and W2 branches are identical for the 60 M|$_\odot$| models since their mass-loss rates are above the de-coupling limit (see Section 5.1 above). In the W2 scheme, all models at |$Z = 10^{-5}$| and |$Z = 10^{-6}$| reach critical rotation. Thus, we find that critical rotation becomes inevitable below |$Z = 10^{-4}$| if we assume zero mass-loss rates when the de-coupling limit is reached. At this point of the evolution, the models generally experience convergence problems, and their further evolution strongly depends on numerical choices. For this reason we do not discuss the model behaviour after reaching critical rotation. In principle, stronger magnetic fields work to prevent high rotation. However, we assume that some non-zero mass-loss rate is required for this process to operate. The magnetic field dependence on the maximum angular velocities shows clearly in the W1 case. While strong magnetic fields are very efficient to brake the surface rotation at high-metallicity, when lowering the metallicity, they become less efficient – even if some (low) mass-loss rates are applied in the models.

Meynet & Maeder (2002b) also studied |$M_{\rm ini} =60$| M|$_\odot$| models at |$Z = 10^{-5}$|⁠. The main difference in their model is the advective treatment of angular momentum transport, which, without including magnetic fields, results in growing radial differential rotation. Although, this is qualitatively different, it still aids outward angular momentum transport in the model (Maeder & Meynet 2001). Combined with the lower mass-loss rates for lower metallicity (scaled down from solar-metallicity conditions, thus similar to our W1 scheme), their model also predicts a spin up early on in the main sequence and reaches critical rotation close to the end of the main sequence (their figs 5 and 7). In our case, we mitigate reaching exactly the critical velocity. Overall, our results are in very good agreement.

6 DISCUSSION

6.1 Ejected mass of chemical species during the main sequence

Understanding chemical enrichment on local, galactic, and even cosmological scales requires well-constrained stellar models. However, the different assumptions lead to large deviations in model predictions.

We demonstrated in Fig. 5 that the He and CNO abundances can be an order of magnitude different when using various mixing efficiencies and magnetic field strengths. Here, we calculate the ejected mass of chemical species over the main-sequence evolution from our models, following the parametrization of Hirschi, Meynet & Maeder (2005) and Higgins et al. (2023).

where |$X_{\rm surf}$| is the surface abundance (mass fraction) of a given isotope, and ZAMS and TAMS are the zero age and terminal age main sequence, respectively. Our results show that the ejected mass is the most sensitive to metallicity and magnetic field strength (Fig. 11). In general, the more efficient Mix2 chemical mixing scheme results in higher nitrogen yields than the Mix1 scheme; however, the carbon and oxygen yields show more complicated trends with chemical mixing scheme. Assuming the W1 wind scheme, the decrease in metallicity results in a gradual decrease in the ejected mass. At |$Z=10^{-3}$|⁠, the strongly magnetized models predict the least amount of ejected mass. On the one hand, these models have less efficient mixing due to their slower rotation. This keeps the surface abundances of carbon and nitrogen relatively unchanged, whereas other models allow for carbon depletion and nitrogen enrichment following the changes in the stellar core. The ejected mass of oxygen is also the smallest for the most strongly magnetized model in a given scheme at a given metallicity. On the other hand, these models also predict shorter main-sequence lifetimes. Since the yields are integrated over a time period, this also means lower values of |$M_{e}$| for the magnetic models. At |$Z=10^{-6}$|⁠, these trends are mostly maintained. In the INT scheme, a very strong, 50 kG magnetic field could have a drastic impact of several orders of magnitude on the ejected mass. The models with different initial magnetic fields strengths in the SURF scheme maintain smaller differences.

In previous works, the rotational properties of low-metallicity massive stars were considered to study the impact on CNO yields by Meynet & Maeder (2002a, b) and Meynet, Ekström & Maeder (2006). We have now demonstrated that additional physical ingredients in the model calculations are able to change the quantitative results already on the main sequence. Interestingly, even if the star initially rotates fast, the nitrogen abundance is not necessarily enhanced when the magnetic field is strong. In particular, the high nitrogen abundance observed in GN-z11 is also accompanied by a relatively low oxygen abundance, which suggests that the source is not a supernova ejecta (which contains oxygen). We note however that, in general, supernova ejecta also contributes to the chemical evolution of the Universe in addition to stellar winds (e.g. Chiaki et al. 2018). While we only focus on main-sequence models in this work, we discuss some potential implications for later evolutionary stages in the next section.

6.2 Implications for low-metallicity supernovae

The properties of the stellar core are important indicators of the star’s subsequent evolution. We evaluate the convective core masses in our models at the TAMS to establish its parameter-space dependence. At this stage, the composition of the core is essentially only helium.

Fig. 12 shows our results combined across various initial magnetic field strengths, the metallicity ranges, the angular momentum transport and magnetic braking schemes, and the two chemical mixing schemes. In initially 60 M|$_\odot$| models, the convective core mass at the ZAMS is |$\approx$|43 M|$_\odot$|⁠. Stars following a classical evolution in the HRD have a convective core that decreases during the main sequence. Overshooting and rotational mixing are able to extend the core size. Since overshooting is assumed with a fixed value in our models (c.f. Paper IV), the main difference is the strength of rotational mixing in the radiative envelope. To a large extent, this is controlled by the Mix1-Mix2 chemical mixing schemes, and the rotational evolution of the star, constrained by magnetic braking. In the limit of QCHE, the core mass can increase during most of the main-sequence evolution.

We find the following main trends. Models using the INT braking scheme tend to have smaller core masses compared to models using the SURF braking scheme. This is directly the consequence of the INT models braking the rotation of the entire star. Hence, they rotate slower in the layers near the core, which makes the core-boundary mixing less efficient. Models utilizing the Mix1 scheme are less efficient in chemical mixing. Therefore, the predicted core masses in the Mix1 scheme are fairly similar. The Mix2 models tend to have much larger helium cores than Mix1 models when the magnetic field is weak and the metallicity is high. Indeed, in general, weaker magnetic fields and lower metallicity also lead to higher core masses. Interestingly, in the Mix1 chemical mixing scheme (both for INT and SURF magnetic braking scenarios), the highest core mass is reached at |$Z=10^{-4}$| (cross symbols). This is unlike the Mix2 case, where the core mass increases as a function of lowering metallicity.

The most striking result of this parameter study is the quantitative difference between the models. For example, for a given initial magnetic field strength or initial metallicity, the helium core masses can differ by up to about 30 M|$_\odot$| for an initially 60 M|$_\odot$| stellar model when considering various modelling assumptions.

The core mass typically changes moderately during the core helium burning phase. For example, recent models of Laplace et al. (2021) indicate that binary effects might slightly increase the helium core mass, while in typical single-stars a slight decrease is expected. Although we do not follow the post-main sequence evolution, the TAMS helium core mass may give some hints towards the expected CO-core masses. In general, the explodability of these models as typical core-collapse events is difficult to predict and will require future investigations until the pre-supernova stages. The onset of pair-instability supernova (PISN), which requires high CO-core masses⁵, may be favoured in low-metallicity environments (e.g. Langer 2012). We find that strongly magnetized models tend to have smaller He core masses (and thus likely smaller CO core masses). Therefore, in addition to low-metallicity environments, a weak magnetic field would also favour reaching the PISN regime in the parameter space considered in this work.

7 FUTURE WORK

Extending high-metallicity stellar evolution model calculations to low-metallicity environments is not straightforward. We have demonstrated that changes in the wind, rotational, and magnetic properties largely affect the model predictions. Thus far the physics of stellar winds and magnetospheres are mostly constrained in the local Universe. If line-driving becomes inefficient to launch the wind, unusual conditions may arise. A low wind density implies a high magnetic confinement for a given strength of a magnetic field. To some extent, the star then becomes a conservative system in terms of mass and angular momentum. For this reason, magnetic and non-magnetic model predictions become identical under the assumption of the W2 wind scheme (e.g. the lowest metallicity models in Fig. 10). If the decrease in wind strength and magnetic braking are gradual, then we expect systematic shifts in the physical quantities, such as surface abundances, ejected masses, and rotational velocities. In the future, we plan to apply sensible considerations and extend this work to zero-metallicity Population III stars.

In order to help constrain primary nitrogen enrichment and study the pre-supernova stages with implication to low-metallicity supernovae, gamma-ray bursts, etc., computations until advanced burning stages are required. Amongst several challenging factors, the stability and evolution of fossil magnetic fields remain elusive beyond the main sequence. A point of interest in this context is the study of convective expulsion (Gough & Tayler 1966; Spruit 1979; Schüssler & Knölker 2001). Namely, if strong convection develops in a stellar layer, the magnetic field lines cannot reconnect, leading to the removal of magnetic flux from a given region. Recently, such stability criterion was studied in one-dimensional stellar evolution codes (MacDonald & Petit 2019; Jermyn & Cantiello 2020). Since extended convective layers develop in the stellar envelope during the post-main sequence, it is an important to consider this when calculating models until advanced burning stages.

8 CONCLUSIONS

In this work, we extend our formalism to calculate surface fossil magnetic field effects in massive star evolutionary models at metallicities lower than that of the SMC. We explore a relevant parameter space to map out key uncertainties in evolutionary model predictions. Our main conclusions are as follows:

• The separation between classical and quasi-chemically homogeneous evolution, as well as the evolution of surface abundances strongly depend on metallicity, the chemical mixing scheme, and the magnetic field strength. Stronger mixing, lower metallicity, and weaker magnetic fields favour the development of QCHE.

• Generally, the surface magnetic field strength weakens by an order of magnitude on the main sequence, following the increase of the stellar radius. However, during QCHE, the stellar radius can remain almost unchanged, leading to a surface magnetic field strength that is constant in time.

• In diffusive angular momentum schemes, the speed of rotation is the driver of chemical mixing in radiative regions. Similar to high-metallicity environments, the rotational evolution is drastically impacted by magnetic braking. However, following our assumptions, this process becomes weaker as a function of metallicity. Instead, the QCHE channel might lead to the onset of Wolf–Rayet type winds, which are able to brake the surface rotation by hundreds of km s|$^{-1}$| on the main sequence.

• We systematically calculate the critical de-coupling mass-loss rates in our models. This shows that already at metallicities of |${\it Z}=10^{-4}$| and |$Z=10^{-5}$| the models with initially 20 and 60 M|$_\odot$|⁠, respectively, may not be able to launch a wind. If we were to assume that there is no mass-loss in such a case (W2 scheme), then the models would tend towards critical rotation.

• The ejected mass of chemical species shows a significant dependence on magnetic field strength. Generally, the stronger the magnetic field, the lower the amount of ejected mass. A strong magnetic field slows the rotation and thus suppresses the surface nitrogen enrichment even if the star initially rotates fast or mixes material efficiently.

• The core masses by the end of the main sequence can vastly differ, depending on the chemical and rotational evolution. For an initial mass of 60 M|$_\odot$|⁠, we find helium core masses in the range of 22–52 M|$_\odot$|⁠. We expect that in low-metallicity environments pair-instability supernova are favoured from initially non- or weakly magnetized progenitors.

Our study is an exploration of physical ingredients in massive star evolutionary models that are unconstrained at low metallicities. As such, we demonstrate that the interpretation based on stellar models is uncertain and requires careful consideration of the underlying assumptions. For this reason, the interpretation of nitrogen excess of distant galaxies could strongly depend on specific stellar progenitor models. Our results for the ejected mass of nitrogen from main-sequence single-star models show a large scatter of several orders of magnitude, highlighting the need to better understand physical processes in massive stars at low metallicities.

ACKNOWLEDGEMENTS

We thank the anonymous referee for providing us with constructive comments that led to improving the manuscript. We thank the mesa developers for making their code publicly available. ZK acknowledges support from JSPS Kakenhi Grant-in-Aid for Scientific Research (23K19071). GC and HN are supported by the NINS International Research Exchange Support Program. HN is also supported by Grant-in-Aid for Scientific Research (23K03468). Au-D acknowledges NASA ATP grant number 80NSSC22K0628 and support by NASA through Chandra Award number TM4-25001A issued by the Chandra X-ray Observatory 27 Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-03060. Numerical computations were carried out on the PC cluster at the Center for Computational Astrophysics, National Astronomical Observatory of Japan.

DATA AVAILABILITY

A full reproduction package is available on Zenodo at https://doi.org/10.5281/zenodo.13309019.

Footnotes

REFERENCES

APPENDIX: mesa MICROPHYSICS

The mesa EOS is a blend of the OPAL (Rogers & Nayfonov 2002), SCVH (Saumon, Chabrier & van Horn 1995), FreeEOS (Irwin 2004), HELM (Timmes & Swesty 2000), PC (Potekhin & Chabrier 2010), and Skye (Jermyn et al. 2021) EOSes.

Radiative opacities are primarily from OPAL (Iglesias & Rogers 1993, 1996), with low-temperature data from Ferguson et al. (2005) and the high-temperature, Compton-scattering dominated regime by Poutanen (2017). Electron conduction opacities are from Cassisi et al. (2007).

Nuclear reaction rates are from JINA REACLIB (Cyburt et al. 2010), NACRE (Angulo et al. 1999), and additional tabulated weak reaction rates (Fuller, Fowler & Newman 1985; Oda et al. 1994; Langanke & Martínez-Pinedo 2000). Screening is included via the prescription of Chugunov, Dewitt & Yakovlev (2007). Thermal neutrino loss rates are from Itoh et al. (1996).