The effects of surface fossil magnetic fields on massive star evolution: IV Grids of models at Solar, LMC, and SMC metallicities

mass-loss

The observed surface magnetic fields of hot, massive and intermediate-mass stars do not show any apparent correlation with stellar and rotational parameters unlike in lower-mass (<2 M ), cool stars (<10 kK), where magnetism due to surface convection and differential rotation ubiquitously produce dynamo activity (Donati & Landstreet 2009;Neiner et al. 2015). Consequently, the organised, large-scale magnetic fields of hot stars are expected to be of fossil origin (Cowling 1945;Spitzer 1958;Mestel 1967Mestel , 2003Moss 2003;Mestel & Moss 2010;Ferrario et al. 2015;Braithwaite & Spruit 2004. The exact origin of observed magnetic fields remains debated. In about 10% of intermediate-mass Herbig Be/Ae stars (which 10% is thought to be the precursors of main sequence Bp/Ap stars), large-scale surface magnetic fields are already observed on the pre-main sequence (Stȩpień 2000;Alecian et al. 2009Alecian et al. , 2013Villebrun et al. 2019;Lavail et al. 2020), which may be acquired from the star-forming disk or generated via a dynamo action inside the star during a fully convective pre-main sequence phase (e.g. Moss 2003;Braithwaite 2012). In addition to magnetic fields possibly remaining from the star formation or pre-main sequence phases, stellar mergers could also amplify seed magnetic fields to a strength sufficient to be detectable (Ferrario et al. 2009;Wickramasinghe et al. 2014;Schneider et al. 2016Schneider et al. , 2019Schneider et al. , 2020, suggesting that there may exist multiple channels to generate globally-organised, 1 However, not all magnetic O-type stars belong to this class (Donati et al. 2002;Petit et al. 2013Petit et al. , 2017Grunhut et al. 2017).
The nature of fossil fields is fundamentally different from contemporaneously generated dynamo fields by a mechanical source (such as convection or differential rotation). Fossil field evolution is purely dissipative with no active field generation counteracting its slow dissipation (Braithwaite & Spruit 2017). In stellar layers where large-scale fossil fields spread through, it is expected that solid-body rotation will develop (e.g. Mestel 1999). In those stellar layers 2 the mechanical source of differential rotation is absent, and consequently small-scale dynamo fields in radiative stellar layers cannot be induced (Spruit 2004). The Tayler instability (Tayler 1973;Goldstein et al. 2019), for example, cannot develop if the radial rotation profile is completely flat, which means that the Tayler-Spruit (or "Ω-type") dynamo cannot be induced in the presence of a fossil field (e.g., Spruit 2004). In fact, while this type of dynamo mechanism in radiative stellar layers was proposed by Spruit (2002), there remains ongoing debate about the necessary electromotive force to operate the dynamo cycle . The simulations of Zahn et al. (2007) suggest that this dynamo cycle does not operate. Despite the contradictory numerical results and the lack of direct observational evidence, dynamos in radiative stellar layers are commonly accounted for in evolutionary models (e.g., Spruit 2002;Maeder & Meynet 2003Maeder 2009;Heger et al. 2005;Yoon et al. 2006;Denissenkov & Pinsonneault 2007;Potter et al. 2012b;Quentin & Tout 2018;Fuller & Ma 2019;Takahashi & Langer 2021). We emphasise that these implementations are not suitable (at least directly) to model stars that are known to host fossil fields.
The time evolution of fossil magnetic fields also remains an unresolved problem. Observed samples of magnetic A-type stars and compact remnants are consistent with the magnetic flux being conserved over time (e.g., Landstreet et al. 2007Landstreet et al. , 2008Wickramasinghe & Ferrario 2005;Neiner et al. 2017;Martin et al. 2018;Sikora et al. 2019a), whereas other observational evidence (including that for OB stars) suggests magnetic flux decay (e.g., Fossati et al. 2016;Shultz et al. 2019b). Fossil magnetic fields are expected to evolve only by Ohmic dissipation (Wright 1969;Spruit 2004;Braithwaite & Spruit 2017), which has a longer timescale than the main sequence nuclear timescale (Cowling 1945;Spitzer 1958). However, Ohmic dissipation remains riddled with uncertainties depending on the exact value of magnetic diffusivity (e.g., Charbonneau & MacGregor 2001) and the geometry of the magnetic field since more complex fields dissipate faster.
Despite the uncertainties regarding the origin and evolution of fossil magnetic fields, it is now well established that they lead to various changes in stellar structure and evolution (e.g., Mestel 1989Mestel , 1999MacDonald & Petit 2019;Jermyn & Cantiello 2020). Two main surface effects, mass-loss quenching and magnetic braking (discussed in detail below), have been shown to drastically modify evolutionary model predictions (e.g., Meynet et al. 2011;Keszthelyi et al. 2017aKeszthelyi et al. , 2019Keszthelyi et al. , 2020Keszthelyi et al. , 2021. For instance, heavy stellar-mass black holes and pair-instability supernovae could be formed from magnetic progenitors even at solar metallicity Georgy et al. 2017).
Thus far, surface magnetic fields have only been detected in Galactic stars. Currently, high-resolution spectropolarimeters used for stellar magnetometry are employed on 4m-class telescopes, which limits observations to bright nearby stars. Using lowresolution spectropolarimetry, Bagnulo et al. (2017Bagnulo et al. ( , 2020 searched for strong magnetic fields in the Magellanic Clouds through the Zeeman effect, which did not lead to definite detections in any of the targets. While high-resolution spectropolarimetry remains largely limited to a Galactic environment, very extensive spectroscopic campaigns in the Magellanic Clouds -in addition to the identified Of?p stars (e.g., Walborn et al. 2015;Munoz et al. 2020) -suggest that the nature of some stars may be explained by invoking surface magnetic fields (Hunter et al. 2008;Brott et al. 2011a;Rivero González et al. 2012;Potter et al. 2012b;Grin et al. 2017;Dufton et al. 2018Dufton et al. , 2020Ramachandran et al. 2021). Observations of known magnetic stars in the Galaxy are often compared to evolutionary models that do not include fossil field effects (e.g., those of Brott et al. 2011a;Ekström et al. 2012;Chieffi & Limongi 2013;Choi et al. 2016;Costa et al. 2021;Grasha et al. 2021), possibly making inferences of stellar parameters rather uncertain. In turn, this can largely impact the derived ages of individual stars and bias isochrone fitting of stellar clusters. Therefore, there is a need for stellar evolution models (and model grids), which take into account mass-loss quenching and magnetic braking (although see Potter et al. 2012b for the latter), thereby affecting detailed evolutionary model predictions and population synthesis studies in both Galactic and extragalactic environments.
The motivation of this study is to help to resolve these issues by presenting and studying an extensive grid of stellar structure and evolution models with metallicities typical of environments in the Solar neighbourhood 3 , Large Magellanic Cloud (LMC), and Small Magellanic Cloud (SMC) that include the effects of surface fossil magnetic fields. The model computations are open source and the entire library of models is available to the community via Zenodo at https://doi.org/10.5281/zenodo.7069766.
This paper is part of a series in which we aim to explore the effects of surface fossil magnetic fields on massive star evolution. In the first paper of the series (Keszthelyi et al. 2019, hereafter Paper I), we used the Geneva stellar evolution code (Eggenberger et al. 2008;Ekström et al. 2012;Georgy et al. 2013;Groh et al. 2019;Murphy et al. 2021) and discussed the mutual impact of magnetic mass-loss quenching, magnetic braking, and field evolution on a typical massive star of initially 15 M at solar metallicity. We studied both solid-body and differentially-rotating models and evaluated their key evolutionary characteristics, showing that strong surface nitrogen enrichment is expected for magnetic models with differential rotation. In the second paper (Keszthelyi et al. 2020, hereafter Paper II), we elaborated on the implementation of massive star magnetic braking in the software instrument (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019 and detailed the magnetic and rotational evolution of the models by performing a parameter test in initial mass, magnetic field strength, and rotational velocity space with 35 models. Then, 72 tailored models were compared with a sample of observed magnetic B-type stars from Shultz et al. (2018Shultz et al. ( , 2019a. A key finding of Paper II is that magnetic stars 3 Since we follow the elemental abundance determinations given by Przybilla et al. (2008Przybilla et al. ( , 2013, Nieva & Przybilla (2012) and Asplund et al. (2009), we refer to this set of models as Solar and not Galactic. Schematic illustration of the magnetic field geometries (solid lines) considered in this study. These geometries are not directly implemented in our 1D models, instead they are used to constrain the corresponding scaling relations. Colours indicate the logarithmic angular velocity; the rotation axis is shown with dotted lines. Left: INT scheme, representative of a dipolar field configuration leading to internal magnetic braking, spinning down all layers of the star. As a phenomenological picture, we assume that convective expulsion (see Appendix C) would not allow the field to relax in the stellar core; however, we assume that the braking can propagate uniformly (including the core). The angular rotation is uniform throughout the star, enforced by a high diffusivity attributed to the dipole field. Right: SURF scheme, representative of a more complex field geometry leading to surface magnetic braking. The magnetic energy can be stored in higher-order spherical harmonics, and twisted field lines only penetrate to some extent of the envelope which we assume to be 20 per cent of the mass fraction. In those layers efficient angular momentum transport is present but we assume that core-envelope coupling is not achieved. Magnetic braking is only applied to those upper layers (see Section 2.6.1). To be able to incorporate the effects of such fields into 1D models, we assume that a quadrupole scaling may be used for calculating the Alfvén radius. could originate from ZAMS progenitors with a variety of parameter combinations. In Paper III (Keszthelyi et al. 2021) we focused on the scenario that some magnetic stars may originate from rapidlyrotating progenitors at the ZAMS, and specifically applied it to the case of the magnetic early B-type star Sco. We found that for this star the simultaneous nitrogen enrichment and slow rotation poses a significant challenge for single-star evolution.
The paper is organised as follows. In Section 2, we detail the assumptions and input physics used in the models. In Section 3, we present and scrutinise the stellar structure and evolution models from our computations. In Section 4, we discuss implications and future work. Finally, we conclude our findings in Section 5.

General strategy
In this work, we follow the general strategy of adopting suitable parametric prescriptions to model the effects of fossil magnetic fields, similar to the approaches presented in Paper I, Paper II, Paper III and references therein. While one-dimensional magnetohydrodynamic (MHD) approaches are possible, they have mostly been developed for dynamo models (e.g., Potter et al. 2012b;Feiden & Chaboyer 2012Quentin & Tout 2018;Takahashi & Langer 2021) and as such are not directly applicable to model fossil fields of intermediate-mass and massive stars. In particular, magnetic transport equations have been developed and used previously in the context of dynamo-generated fields (e.g., Spruit 2002;Maeder & Meynet 2003Heger et al. 2005;Yoon et al. 2006;Potter et al. 2012b;Kissin & Thompson 2018;Quentin & Tout 2018;Fuller & Ma 2019;Takahashi & Langer 2021). Although the characteristics (for example, the scale and time evolution) of dynamo models are incompatible with those of fossil fields (see Section 1), the transport equations may follow similar implementations. Here, we opt to artificially increase the diffusivity instead of testing "magnetic" transport equations (Section 2.6.2). Such equations would introduce more free parameters and further assumptions regarding the geometry, structure, and radial dependence of the magnetic field. Clearly, further research and observational verification is required before an appropriate one-dimensional magnetic transport process could be reliably incorporated to model stellar evolution with fossil magnetic fields (however, see Schneider et al. 2020).  and  presented a comprehensive approach applicable for fossil fields, showing however that the impact on hydrodynamic equilibrium and energy transport are modest even for strong magnetic fields. To this extent, it is indeed appropriate to use parametric prescriptions and focus on the major, measurable effects 4 that fossil magnetic fields have, namely changing the mass loss (Section 2.5) and rotation (Section 2.6) of the star, affecting chemical mixing and angular momentum transport.
One of the major modelling challenges is that the geometry and alignment of the magnetic field play a significant role in the corresponding physical description. It has been demonstrated that a seed magnetic field can relax into a stable axisymmetric (around the magnetic axis) configuration if the magnetic flux is centrally concentrated, or into a non-axisymmetric (around the magnetic axis) configuration otherwise (Braithwaite & Spruit 2004;Braithwaite & Nordlund 2006;Braithwaite 2008). In both cases, the latitudinal averaging is inappropriate to model the magnetic field in 1D 5 .
For simplicity, we assume that the field is aligned with the rotation axis of the star since appropriate scaling relations for oblique fields (tilted with respect to the rotation axis) are still in development. However, the obliquity angles inferred from observations appear to follow a random distribution, which suggests that, apart from a few possible exceptions, massive stars generally possess magnetic fields that are inclined with respect to the rotation axis (e.g., Khalack et al. 2003;Shultz et al. 2019a;Sikora et al. 2019a). Recent work from ud-Doula (2020) suggests that oblique rotation leads to decreasing the efficiency of magnetic braking. This effect could be incorporated in our models via a suitable scaling factor in future studies. However, the efficiency of magnetic braking, in the evolutionary context, is also largely dependent on magnetic field evolution, which still needs to be better constrained (e.g., Paper III). 4 See e.g., Driessen et al. (2019a) for mass-loss quenching, and e.g. Townsend et al. (2010); Oksala et al. (2012); Song et al. (2021) for magnetic braking. 5 For example, in the case of an axisymmetric dipole geometry, both poloidal and toroidal components must exist in the stellar interior (Braithwaite & Spruit 2004). The poloidal field strength and orientation relative to the normal to the surface varies over latitudes. The toroidal field, confined by closed poloidal lines, has zero strength along the polar rotation axis and reaches its maximum along the equatorial plane (see, e.g., Figure 4 of Braithwaite 2008). A latitudinal averaging instead assigns a mean value to the poloidal and toroidal field components along a radius.
In fact, magnetic field evolution is closely tied to the question of the field geometry and, in this regard, new insights are gained from extensive monitoring campaigns, which can reveal the surface properties of magnetic fields. Although a purely dipolar field geometry generally matches observations (Grunhut et al. 2017;Shultz et al. 2018), modest deviations from pure dipolar geometries are now identified (Leto et al. 2018;Das et al. 2020;David-Uraz et al. 2021a). In other cases, contributions from higher-order harmonics are also identified (e.g., Shultz et al. 2018;Kochukhov et al. 2019); however, the dipole is the strongest component, which consequently drives the main physical effects. In a few cases, observations have also identified stars with uniquely complex magnetic fields, which cannot be described with a dominant dipole component (e.g., Sco, Donati et al. 2006;Kochukhov & Wade 2016;Shultz et al. 2018). In these cases, most of the magnetic energy is stored in higher-order spherical harmonics, although a weak dipole contribution may still be present.
Since the structure of large-scale magnetic fields in stellar interiors is still considerably uncertain, we aim to test two limiting cases -illustrated in Figure 1 -to apply depth-dependent magnetic braking (c.f. Paper II). In the models with internal magnetic braking (INT, described in detail below), we assume that the magnetic field is dipolar. We further assume that the field is present in the entire stellar envelope and is able to achieve core-envelope coupling (however, only as a phenomenological picture, we assume that it is excluded from the stellar core). In nature, the stellar envelope must also have a closed toroidal field for a stable configuration (Wright 1969(Wright , 1973Braithwaite & Spruit 2004;Braithwaite & Nordlund 2006;Akgün et al. 2013), which we cannot directly take into account in our 1D models. The magnetic field is curl-free above the stellar surface; any toroidal field diffuses to a poloidal structure. We also introduce a set of models with surface magnetic braking (SURF, see also below) as a limiting scenario to contrast the INT models with. The complex magnetic field geometries clearly cannot be translated to our 1D models, therefore we make various simplifying assumptions (discussed in detail throughout Section 2). We reiterate that the field geometry cannot be directly included in our models, only via the scaling relations. The limitations related to the 1D parametrisation can likely only be resolved once 2D stellar evolution modelling becomes feasible (Espinosa Lara & Rieutord 2011;Lovekin 2020;Reese et al. 2021).
Angular momentum is always lost from the outermost layers of the star. We use the INT/SURF schemes to control the propagation of this loss to the stellar interiors along with efficient angular momentum transport. In the INT models, all stellar layers efficiently transport and lose specific angular momentum. Motivated by the simulations of Braithwaite (2008), we assume that in the SURF models the magnetic-field driven angular momentum transport and rotational braking only affect the upper 20% of stellar mass. Although even for complex surface magnetic fields there may be weak dipole contributions, we neglect here any possible magnetic angular momentum transport in the deep stellar interiors so that we are able to test a limiting case in which radial differential rotation may develop in regions of the star where the magnetic flux is assumed to be negligible. Furthermore, to generalise from the range of possible field configurations deduced from observations, we assume that at least the Alfvén radius (see section 2.4) has to be smaller in the SURF models than in the INT models for a given surface field strength. For this reason, we use a quadrupole scaling in the SURF models to obtain the Alfvén radius. This results in less efficient magnetic braking for complex fields compared to dipole-dominated geometries. The field geometry defines the wind flow and mass flux from the stellar surface. We evaluate in Appendix B how an actual quadrupole field geometry would affect mass-loss quenching; however, since this effect only concerns the highest-mass models, for simplicity we adopt a formalism where only the Alfvén radius is changed in the SURF models compared to the INT models (see Section 2.5).

General model setup
We use the software instrument Modules for Experiments in Stellar Astrophysics release 15140 (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019 and Software Development Kit (SDK) version 20.12.1 , and carry out our computations on the Dutch supercomputers Cartesius 6 and Snellius 7 . In this study, we consider main sequence models with initial masses from 3 to 60 M (Table 1). This choice is made since there are still considerable uncertainties in stellar evolution, as well as in magnetic field evolution on the post-main sequence. The computations begin by relaxing the initial model and we consider the Zero Age Main Sequence (ZAMS) to begin at the time when the initial abundance of core hydrogen has decreased by 0.3 percent. The endpoint of the models is the Terminal Age Main Sequence (TAMS), which we consider at the time when the core hydrogen mass fraction drops below 10 −5 .
The microphysics are summarised in Appendix A. To calculate the nuclear reaction rates, we use 's default "basic.net" option. To model convective mixing, we assume a mixing length parameter of MLT = 1.8 (e.g., Canuto & Mazzitelli 1992;Canuto et al. 1996) and the Henyey formalism (Henyey et al. 1965) in . Common values of MLT range from 1.5 to 2.0, mostly based on solar and solar-type star calibrations (e.g., Bonaca et al. 2012). For simplicity, we assume the same value in all models, although some studies predict mass dependence (e.g., Yıldız et al. 2006). Semiconvective and thermohaline mixing are not used 8 . Convective boundaries are determined via the Ledoux criterion. Core overshooting is applied in the exponential scheme with parameters ov = 0.015 and 0 = 0.005, which roughly corresponds to a step overshoot parameter ov of 10 percent of the local pressure scale height , see further discussion in Appendix C). ov may be mass dependent (Castro et al. 2014). Commonly used values of ov range from 0.1 to 0.335, (Schaller et al. 1992;Brott et al. 2011a). Exponential overshooting at non-burning convective regions is adopted with ov = 0.0010 and 0 = 0.0005. 's MLT++ scheme is not applied (see e.g., Poniatowski et al. 2021).
We employ high spatial and temporal resolution by setting mesh_delta_coeff = 1. and time_delta_coeff = 1., in addition to setting varcontrol_target = 1.d-4. This results in an average of 2000-3000 zones for our stellar structure models, and the evolutionary models consist of a hundred to a few thousand structure models (each corresponding to one time step), mostly depending on the initial mass. 6 https://userinfo.surfsara.nl/systems/cartesius 7 https://userinfo.surfsara.nl/systems/snellius 8 Neither of these mixing mechanisms is expected to significantly change main sequence models. It was shown by Charbonnel & Zahn (2007) that a strong magnetic field could inhibit thermohaline mixing in descendants of Ap stars (see also, e.g., Denissenkov et al. 2009). On the other hand, Harrington & Garaud (2019) finds that thermohaline mixing could be enhanced by an aligned magnetic field.
We consider one initial rotation rate 9 in all models by relaxing an initial ratio of Ω/Ω crit = 0.5. In our models with solar metallicity, this approximately corresponds to an initial equatorial surface rotational velocity (surf_avg_vrot in ) of 300 to 370 km s −1 in the initial stellar mass range from 3 to 60 M . The critical angular velocity is adopted as defined in star_utils.f90: where the Eddington parameter is Γ = rad / Edd , G is the gravitational constant, and and are the mass and radius taken at the photosphere. This definition only plays a role in setting the initial rotational velocity 10 . The initial model relaxes from solid-body rotation to a new configuration constrained by angular momentum distribution and transport.
In Paper II we found that a lower rotation rate (initially Ω/Ω crit = 0.2) leads to smaller differences between models with and without magnetic fields. This is simply because magnetic braking is less efficient when the rotation is slow 11 . In Paper II, we showed that at least a few magnetic stars were best matched with models that had an initial rotation rate of Ω/Ω crit = 0.8. If the initial rotation rate was higher than considered in this study, it would i) alter the early evolution of the models on the Hertzsprung-Russell diagram (HRD, as shown in Figure 3 of Paper II), and ii) impact the quantitative predictions regarding the surface chemical enrichment. It is worth noting that the most rapidly rotating (presumably young) magnetic stars have present-day surface rotational velocities of about 300 km s −1 , which is close to 50% of the critical rotation defined in our Equation 1 (e.g., Oksala et al. 2010;Grunhut et al. 2012;Shultz et al. 2019b;Song et al. 2021). Nevertheless it remains unknown what initial rotation rates could characterise the entire sample of magnetic massive stars.
Finally, given the supporting observational evidence by some Spectroscopic studies have focused on large samples to obtain the distribution of projected rotational velocities in the Galaxy (Howarth et al. 1997;Huang et al. 2010;Simón-Díaz & Herrero 2014;Simon-Diaz et al. 2016;Simón-Díaz et al. 2017;Holgado et al. 2022) and in the Magellanic Clouds (Martayan et al. 2006(Martayan et al. , 2007Ramírez-Agudelo et al. 2013Dufton et al. 2019Dufton et al. , 2020. The findings indicate a Gaussian distribution of sin , with different peak values depending on physical (spectral type, mass, etc.) and observational characteristics (sample size, magnitude limit, etc.). The typical peaks of the distributions are around 100 km s −1 . Considering that this value needs to be corrected for the (usually unknown) inclination angles and that it reflects on the current rotation after a given star or population have evolved away from the ZAMS, it is generally assumed that the canonical initial rotational velocities of massive stars are of the order of 300 km s −1 .
Although for Galactic O-stars the IACOB survey (Simón-Díaz & Herrero 2014;Simon-Diaz et al. 2016;Simón-Díaz et al. 2017;Holgado et al. 2022) has shown somewhat lower values than previous studies (≈ 100-200 km s −1 ), which could be consistent with lower initial rotational velocities. Our chosen input parameter for the initial rotation rate of Ω/Ω crit = 0.5 reflects closely on the canonical value around 300 km s −1 , identified in the sample of Howarth et al. (1997). 10 In fact, in the Eddington luminosity is calculated from the total opacity. Since the precise definition of Ω crit is significantly more complex (see Puls et al. 2008 for further discussion), we stress here that Ω crit is only used as an input option to set the rotational velocity. 11 In Paper I, we also demonstrated that for a typical non-rotating 15 M model at solar metallicity, mass-loss quenching is modest. As shown by Petit et al. (2017), who computed non-rotating solar metallicity models between 40 and 80 M , the evolutionary impact of magnetic mass-loss quenching becomes significant at higher masses. studies (Kochukhov & Bagnulo 2006;Wickramasinghe & Ferrario 2005;Neiner et al. 2017;Martin et al. 2018;Sikora et al. 2019a), we assume magnetic flux conservation (Alfvén's theorem, Alfvén 1942) such that the surface magnetic field strength is obtained from: with eq,ini being the initial surface equatorial magnetic field strength (which is assumed in the models at the ZAMS), while ★,ini and ★ are the initial and current stellar radii. For further discussions on magnetic field evolution, we refer the reader to Paper I, Paper II, Paper III and references therein. We adopt a large range of initial equatorial magnetic field strengths (from 0 to 50 kG; Donati & Landstreet e.g., 2009;Shultz et al. e.g., 2019b.) The run_star_extras file (available as part of a full reproduction package on Zenodo at https://doi.org/10.5281/ zenodo.7069766) is used to modify the wind, torque, angular momentum transport, and chemical element transport routines. The parameters varied in the models are summarised in Tables 1, 2, and 3, and are discussed below.

Metallicity
We compute models for 3 different metallicities (the key elements are summarised in Table 2 -the full list of abundances is available via the model files shared on Zenodo) that are representative of the Solar neighbourhood, LMC, and SMC. For the Solar composition, we assume the hydrogen and helium mass fractions along with a metallicity of = 0.014 (Asplund et al. 2009). The metal fractions are adopted following the works of Przybilla et al. (2008Przybilla et al. ( , 2013 and Nieva & Przybilla (2012), which updated some elements compared to the Asplund et al. (2005Asplund et al. ( , 2009) abundances, considering B-type stars in the solar neighbourhood. The baseline values for the Magellanic Clouds are still subject to ongoing investigations (e.g., Dufton et al. 2020;Bouret et al. 2021). For the helium and metal abundances in the LMC and SMC, we adopt the mean values listed in Tables  5 and 6 of Dopita et al. (2019). These mean values are a result of 9 separate investigations using different approaches, which include atmospheric determinations of hot stars, supernova remnants, and H regions. In all 3 metallicities, we use the Lodders (2003) isotopic ratios. The metallicity adopted for the chemical composition is fully consistent with the metallicity used for the opacity tables.

Alfvén radius
The Alfvén radius characterises a critical distance at which the magnetic energy density and the gas kinetic energy density are equal. Alternatively, it can also be cast as the inverse square of the Alfvénic Mach number. Its definition plays an important role in both mass-loss quenching (Equations 7-8) and magnetic braking (Equations 9-10). For a dipolar field configuration, ud-Doula et al.
with ★ the equatorial magnetic confinement parameter, defined as: with eq the equatorial magnetic field strength, =0 the mass-loss rate in absence of a magnetic field, and ∞ the terminal velocity 12 (ud-Doula et al. 2009). Observations typically reconstruct a polar field strength p from the line-of-sight disc-integrated (so called longitudinal) magnetic field strength (Donati & Landstreet 2009). The equatorial field strength is exactly one half of the polar field strength. We use Equation 3 to obtain the Alfvén radius in the INT models, which we assume to be characterised by a predominantly dipolar field configuration ( Figure 1). In the SURF models we assume the field to be more complex, in which case the definition of the Alfvén radius is non trivial. For the sake of simplicity we assume that the Alfvén radius takes the form of a scaling appropriate for a quadrupole field geometry, such that: following the parametrisation in Equation 9 of ud-Doula et al. (2008). This ensures that for a given field strength A is less in the SURF case than in the INT case, leading to less efficient magnetic braking.

Mass-loss schemes and terminal velocities
The models include mass loss. Even though this is modest for the lower-mass stars (and the driving mechanism is not unambiguously identified as for more massive stars), it can impact their rotational evolution given the longer nuclear timescale. For this reason, we apply commonly used mass-loss rates of hot massive stars also to lower mass main sequence stars in our grid. While the highermass stars typically reach the TAMS at eff > 20 kK, we describe here the detailed treatment implemented in our extension for completeness and to aid further studies focusing on complementing this work with post-main sequence models. For massive stars with eff > 10 kK, the mass loss is powered by radiative line driving (e.g., Lucy & Solomon 1970;Castor et al. 1975;Puls et al. 2008). In this regime, we apply the rates derived by Vink et al. (2000Vink et al. ( , 2001, decreased by a factor of 2 for all models in the 3 -60 M range for consistency. The choice to reduce the nominal mass-loss rates is motivated by the growing evidence both from observations, suggesting that mass-loss rates are lower when accounting for wind clumping (e.g., Bouret et al. 2005;Fullerton et al. 2006;Trundle & Lennon 2005;de Almeida et al. 2019;Brands et al. 2022), and from new modelling approaches (e.g., Muĳres et al. 2012;Krtička 2014;Krtička & Kubát 2017;Krtička et al. 2021;Sundqvist et al. 2019;Björklund et al. 2021). When using these rates, we apply metallicity-dependent winds with a scaling of ∼ 0.85 (e.g., Vink et al. 2001;Mokiem et al. 2007).
In agreement with the partitionings in effective temperature, we estimate the terminal wind velocity ∞ via: where , ★ , ★ , and Γ are respectively the gravitational constant, the stellar mass, the stellar radius, and the Eddington parameter for pure electron scattering. The terminal wind velocity is obtained from the escape velocity as a simple step function by adopting ∞ = 2.6, 1.3, 0.7 at eff > 20 kK, 20 kK > eff > 10 kK, and eff < 10 kK, respectively (Lamers et al. 1995;Vink et al. 2000;Kudritzki & Puls 2000). The typical terminal velocities at the ZAMS range from 800 to 3000 km s −1 for models with initial masses from 3 to 60 M , respectively. We calculate the rotational enhancement on the mass-loss rates 13 as described by Maeder & Meynet (2000). This requires defining the difference of the force multiplier parameters (= − , that is, the exponent related to the line-strength distribution function minus the exponent quantifying the change in ionisation balance), which we adopt as a simple step function with values of 0.6, 0.5, 0.4, corresponding to the above-mentioned effective temperature ranges (see Pauldrach et al. 1986;Lamers et al. 1995;Puls et al. 2000). The alternative calculation of rotational enhancement built into is not used (see Paper II Section 3.9 for details). 13 See also the recent study of Brinkman et al. (2021).

Magnetic mass-loss quenching
The overall field configuration that extends into the wind outflow is governed by the competition between the kinetic energy of the wind and the magnetic energy of the field. The ionised stellar wind material is forced to flow along magnetic field lines. However, as the wind kinetic energy density has a shallower decline than the magnetic energy density, the field loops can only confine wind material up to a certain radius. Within closed field loops, material becomes trapped and eventually falls back onto the surface (unless centrifugally supported). To account for the global, time-averaged effect of the magnetosphere, the mass-loss rates are systematically reduced. Following the works of ud-Doula et al. (2008,2009), the mass-loss quenching parameter B is defined as: and where A , K , and c are the Alfvén radius, the Kepler co-rotation radius, and the closure radius in units of the stellar radius, respectively (see Petit et al. 2017, Paper I, Paper II, Paper III, and references therein). The closure radius, defining the distance from the stellar surface to the last closed magnetic loop, is approximated as c ∼ ★ + 0.7 ( A − ★ ), see ud- Doula et al. (2008). is the mass-loss rate that a non-rotating magnetic star would have. is further scaled by the rotational enhancement rot (specified in Section 3.9 of Paper II) such that the effective mass-loss rate is obtained from eff = B · rot · =0 . The magnetic mass-loss quenching parameter (equivalent to the escaping wind fraction 14 ) can take values between 0 and 1, depending on the magnetic field strength. A strong magnetic field (with a strength of tens of kG) may lead to only a few percent of the wind material actually escaping the star Georgy et al. 2017, Paper I). Let us also note that the conditions in the above equations are equivalent to distinguishing between dynamical magnetospheres (if A < K ) and centrifugal magnetospheres (if K < A ), a classification introduced by Petit et al. (2013).
The use of Equation 8 is a refinement compared to previous implementations. For situations when the Alfvén radius is larger than the Kepler co-rotation radius (centrifugal magnetospheres), the magnetosphere is expected to be less efficient at quenching wind mass-loss compared to dynamical magnetospheres (ud-Doula Table 3. Magnetic braking and chemical mixing schemes

INT:
solid-body rotation, braking the rotation of the entire star SURF: differential rotation, braking the rotation of the surface  2008,2009). This is because material injected by the wind into the centrifugal magnetosphere is not returned to the stellar surface by gravity, but is instead ejected away from the star once the critical centrifugal breakout density is exceeded Owocki et al. 2020). This can lead to substantially larger values of B for a rotating as compared to a non-rotating star, however in practice rapid spin-down means that centrifugal magnetospheres are relatively short-lived and the incorporation of this modification to the mass-quenching prescription does not have a strong effect on evolution. It is generally expected that the evolution proceeds from centrifugal to dynamical magnetospheres (Shultz et al. 2019b, Paper I, Paper II).
In this approach, the magnetic mass-loss quenching parameter is an average quantity. MHD simulations of non-rotating magnetospheres predict up-and down-flows of material varying on short dynamical timescales (ud-Doula & Owocki 2002), which can manifest as stochastic variability in magnetospheric emission lines (ud-Doula et al. 2013), however time-averaged models provide a good reproduction of emission line properties (e.g. Sundqvist et al. 2012;Owocki et al. 2016;Erba et al. 2021) and these short-term, stochastic variations can therefore be confidently neglected over evolutionary timescales. In the case of rapid rotators, 2D MHD simulations led to the expectation that breakout events would be similarly stochastic, leading to emptying of the centrifugal magnetosphere and largescale magnetospheric reorganisation (ud-Doula et al. 2006(ud-Doula et al. , 2008. However, no indication of large-scale changes has been observed Shultz et al. 2020), leading Shultz et al. (2020) and Owocki et al. (2020) to infer that breakout events are characterised by small spatial scales and occur more or less continuously, such that the centrifugal magnetosphere is maintained nearly continuously at the breakout density. As a result, it is therefore appropriate to treat magnetospheric mass-drainage via breakout as an effectively continuous process over evolutionary timescales and apply Equations 7 and 8.

Angular momentum transport and loss
Magnetic fields are much more efficient at transporting angular momentum than purely hydrodynamic processes such as meridional currents and shear instabilities (e.g., Mestel 1999;Spruit 1999Spruit , 2002Kulsrud 2005;Braithwaite & Spruit 2017). The rotation of the star leads to Maxwell stresses, which result in losing angular momentum from the star.
In Paper I, models were used, where magnetic braking is adopted as a boundary condition to internal angular momentum transport, directly affecting the uppermost layer of the stellar models. In Paper II, two kinds of models were introduced to account for the uncertainty regarding how deeply fossil magnetic fields are anchored in massive stars. In the INT models, magnetic braking was applied to the entire star, decreasing uniformly the specific angular momentum in all layers. In the SURF models, magnetic braking was set to remove specific angular momentum from a very near-surface reservoir. In Paper III, models were contrasted with a implementation where magnetic braking was applied to most of the stellar envelope. In , two configurations were used to model internal angular momentum transport: one with only hydrodynamic instabilities correctly accounted for via an advecto-diffusive equation (allowing for shears to develop in deeper layers), and one with a purely diffusive equation, in which solid-body rotation was established. In , we relied on the nominal hydrodynamic transport processes since they are used in a purely diffusive assumption, leading to nearly solid-body rotation on the magnetic braking timescale. Here, we make some further refinements and adjustments compared to these approaches, particularly accounting (indirectly) for the field geometry as depicted in Figure 1.

Magnetic braking
Stellar rotation bends and twists magnetic field lines in the azimuthal direction. Magnetic field lines can transport and store angular momentum, and the associated Maxwell stresses are very efficient at transferring angular momentum to the surrounding plasma. Once the angular momentum is imparted from the field to the gas, the wind material carries it away, leading to a spin-down of the star. This process is commonly referred to as (wind) magnetic braking.
In a pioneering series of works, analytical and numerical MHD simulations were developed, confirming that the Weber & Davis (1967) model (see also , Parker 1958;Mestel 1968) leads to an appropriate scaling relation also for massive stars (ud-Doula & Owocki 2002;ud-Doula et al. 2008ud-Doula et al. , 2009Owocki & ud-Doula 2004;Townsend & Owocki 2005). Following the work of ud-Doula et al. (2009), the total -wind and magnetic field induced -loss of angular momentum can be expressed via: with d B /d the rate of angular momentum loss from the system, Ω ★ the surface angular velocity, and the Alfvén radius (defined in Equations 3 and 5). As this equation accounts for the gas and field driven angular momentum loss (ud-Doula et al. 2009), it yields the angular momentum loss resulting purely from mass loss when surf = 0. As specified in Paper II, we have adjusted the angular momentum lost via mass loss to avoid double counting. In Equation 9, the numerical term 2/3 arises from integrating over latitudes. We note that this equation is not applicable when the effective mass-loss rate, as introduced above, is exactly zero (this situation does not happen in our models). In the strong confinement limit, when B → 0, the effective mass-loss rate can become very small. In this case, a strong magnetic braking can still be achieved since the Maxwell stresses driving the angular momentum transport are independent of the plasma flow. As long as there is wind material at a radial distance larger than the last closed magnetic field line, i.e., the star is not surrounded by vacuum, the field can impart angular momentum to the plasma. In Paper II and Paper III, Equation 9 was implemented into via changing the specific angular momentum in given layers of the star, such that a summation over mass yields the total rate of angular momentum loss as defined in Equation 9. It is coded as: where d B /d is the rate of specific angular momentum change (dubbed as "extra_jdot" in ). The negative sign is added to reduce the reservoir (i.e., to account for loss), d B = (d B /d ) · d is the total angular momentum lost per time d , INT/SURF is the angular momentum reservoir of the entire star (INT) or of defined layers in the stellar envelope (SURF; see Figure 1 and discussion below), is the specific angular momentum of a layer (called "j_rot" in ), d is one timestep in the computation 15 , is an index running through all layers, and is the index of the last layer where magnetic braking is applied. Therefore, Equation 10 indicates how to distribute the total angular momentum lost per unit time (d B /d given by Equation 9) in given stellar layers. Taking the sum of the specific angular momentum lost per unit time d B /d with respect to mass, we recover the left-hand side term.
To distribute the total angular momentum lost per unit time, the summation goes over the layers of the entire star in the INT case ( ≈ 3000 zones), whereas in the SURF case it goes from the photosphere to a lower boundary. This boundary is always in the radiative stellar envelope of our models. However, more massive models have larger convective cores, and thus for very massive stars (> 60 M ), this condition may need to be revised as we do not expect the fossil field to be able to penetrate into the convective core. In the SURF models, ≈ 200 zones undergo magnetic braking. Here, we chose the boundary layer where = / ★ = 0.8 (the enclosed mass is 80% of the total mass) since Braithwaite (2008) demonstrated that complex, non-axisymmetric fields (around the magnetic axis) can form if the magnetic flux is initially not centrally concentrated, leading to a stable magnetic field configuration in which twisted magnetic field lines spread throughout the stellar surface layers. In the simulations of Braithwaite (2008), a strong toroidal field (enclosed by poloidal field lines) is present in approximately 20% of the upper mass fraction and this motivates our choice for this parameter.

Angular momentum transport
The main impact of the angular momentum transport equation in stellar interiors is to change the angular velocity profile Ω( ), which is also measurable via modern asteroseismology (see, e.g., Aerts et al. 2019 for a comprehensive review). In , angular momentum transport is modelled in a fully diffusive scheme. Note that this approach inadequately models the meridional currents 16 , which are an advective process by nature.
solves the angular momentum transport equation following Equation 46 of Heger et al. (2000), which is based on the works by Endal & Sofia (1978) and Pinsonneault et al. (1989), that is: where AM is the total diffusion coefficient responsible for angular momentum transport, while , , and are the radius, density, and enclosed mass, respectively, and is the time.
In the non-magnetic models, we assume that AM is constructed as a sum of four diffusion coefficients (resulting from dynamical and secular shear, meridional circulation, and GSF instability), which are the same as used for the Mix1 chemical mixing scheme in Equation 16; however, not scaled by any efficiency parameters for angular momentum transport. For simplicity and a consistent treatment of angular momentum transport, we also use these diffusion coefficients for angular momentum transport when a different chemical mixing scheme is adopted (Mix2, see below). 15 We use a timestep control, specified in Paper II, which prevents the star model from fully exhausting specific angular momentum in any layer. 16 Meridional currents are large-scale flows arising from the thermal imbalance between the polar axis and the equatorial regions in a rotating star.
In the INT models (see also Figure 1), we assume that the magnetic field is capable of establishing radially uniform (solid-body) rotation throughout the entire star. This is representative of an axisymmetric magnetic field that "freezes" rotation along the poloidal field lines following Ferraro's theorem (Ferraro 1937). We model this by using the controls set_uniform_am_nu_non_rot = .true. and uniform_am_nu_non_rot = 1.d16 such that a high diffusivity ( AM = 10 16 cm 2 s −1 ) leads to efficient angular momentum transport and hence solid-body rotation throughout the entire star. Unfortunately, the naming conventions here are somewhat confusing as these controls are applied to the entire star regardless of the convective/radiative nature of given layers. Otherwise "am_nu_non_rot" refers to layers of the star with convective mixing. The precise value of this quantity is not crucial so long as it achieves solid-body rotation. Above a critical value, the diffusivity can saturate, meaning that an already flat Ω profile will remain unchanged if an even higher diffusivity is applied. The "default" value for this control is AM = 10 20 cm 2 s −1 . Such a high diffusivity would mean a diffusion timescale ( D ≈ 2 / AM ) of a few hours, which is physically not justified. The saturation, i.e., solidbody rotation for a given diffusivity may happen for diffusivities > 10 10 cm 2 s −1 , depending on model specifics such as mass and evolutionary stage.
In the SURF models, we distinguish between three regions of the star i) the stellar core, ii) the envelope from = 0.8 to the stellar core, and iii) the envelope above = 0.8 in which magnetic braking is applied (see above). On the main sequence, the cores of massive stars are convective, dominated by strong turbulent mixing. In , this is modelled by a high diffusion coefficient (relying on mixinglength theory) that establishes a constant angular velocity profile, that is, the core is rigidly rotating. In the radiative layers between the stellar core and the boundary of = 0.8, the usual hydrodynamical instabilities (dynamical and secular shear, meridional circulation, GSF instability) transport angular momentum. More directly, the assumption here is that there is no magnetic coupling between the stellar core and the envelope. While even for complex surface fields there may be weak dipole components in the deep stellar layers which may contribute to angular momentum transport, we neglect those here to be able to test a limiting, boundary case, in which differential rotation may develop between the core and the surface. For a consistent comparison, in both Mix1/Mix2 chemical mixing schemes (see below), we apply the same treatment of angular momentum transport in this region.
The fossil magnetic field may relax into a non-axisymmetric configuration, strongly impacting the upper stellar layers (Braithwaite 2008). In these layers (with 20 per cent of the stellar mass in our models), we apply a high diffusion coefficient of AM = 10 16 cm 2 s −1 via the other_am_mixing subroutine to account for the expected effect of the magnetic field.
In both INT and SURF cases, for layers with increased angular momentum transport attributed to the magnetic field, the angular momentum transport equations are of secondary importance in the sense that we expect an appropriate transport equation to result in a flat angular velocity profile, thereby deviating from non-magnetic models. One would also expect that in those layers where the fossil magnetic field is present, hydrodynamical instabilities could not transport angular momentum.
Further guidance regarding the internal rotation profile and magnetic field properties can also be obtained observationally using (magneto-)asteroseismology (see recently Lecoanet et al. 2022). For instance, radial differential rotation was observed in several massive stars using the rotational splitting of gravity modes (e.g. Aerts et al. 2003;Triana et al. 2015). On the other hand, the nearly identical surface and core rotation of red giant stars requires very efficient transport (e.g., Moyano et al. 2022). Using asteroseismic analysis of Kepler data, it has indeed been attributed to magnetic fields (Fuller et al. 2015). Due to possible mode suppression by strong magnetic fields, magnetoasteroseismology remains an elusive target, having been performed for only a few massive stars (e.g. HD 43317 and V2052 Oph; Briquet et al. 2012;Buysschaert et al. 2018). However, the advent of nearly all-sky high-precision space-based photometry can help further this line of inquiry, with large asteroseismic target lists of OB stars already being assembled (e.g. Burssens et al. 2020).

Rotational mixing of chemical elements
Following Pinsonneault et al. (1989), rotational mixing of chemical elements is commonly applied via the diffusion equation in onedimensional stellar evolution models: where is the mass fraction of a given element , is the time, and are the mass coordinate and mean density at a given radius , chem is the sum of individual diffusion coefficients contributing to chemical mixing (see also Salaris & Cassisi 2017), and the last term accounts for nuclear burning.
In this approach, the central question is how to encapsulate inherently three-dimensional physical processes and apply them via a single parameter chem . In this study, we contrast two commonly used approaches. Spectroscopic studies of massive stars often find discrepancies between the observed and predicted surface abundances from rotating stellar evolution models computed with a given scheme of chemical mixing (e.g., Trundle et al. 2004;Martins et al. 2017;Markova et al. 2018). Recent works suggest that such discrepancies may be resolved by including additional processes in the calculations, for example, internal gravity waves (and magnetic fields) lead to a more complex physical interplay between various processes and a variety of mixing profiles (Aerts et al. 2019;Bowman et al. 2020;Michielsen et al. 2021;Pedersen et al. 2021).

Basic thermodynamic quantities
Before introducing the diffusion coefficients, we briefly outline the most important thermodynamic quantities that enter into those equations. The thermal diffusivity is defined as: where is the radiation constant, the speed of light, the mean radiative opacity, the specific heat capacity per unit mass at constant pressure, the temperature, the density, and the pressure. The different ∇-s below denote the adiabatic, radiative, and chemical composition ( ) gradients: where is the entropy, the mean molecular weight, and the gravitational constant. The local luminosity is the rate of energy transported outward through a sphere of radius , and is the enclosed mass. From Equation 4.22 of Maeder & Zahn (1998), the derivatives from the equation of state are:
To be able to compare our results to common model grids, in the Mix1 scheme we adopt the diffusion coefficient applied in Equation 12 as: where the individual diffusion coefficients are described according to Heger et al. (2000). Transport by dynamo mechanisms and by the Solberg-Høiland instability are not considered as their contribution to chemical mixing has been debated (e.g., Yoon et al. 2006;Brott et al. 2011a). Of particular interest is the meridional circulation term, which was described via the circulation velocity in the radial direction by Kippenhahn (1974) and constructed into a diffusion coefficient by Endal & Sofia (1978). However, the base formulation of the problem in terms of a steady-state circulation by Vogt (1925), Eddington (1925) and Sweet (1950) has been disputed by, e.g., Busse (1981Busse ( , 1982; Zahn (1992) -see further discussion by Rieutord (2006). The simple summation of the various processes by Heger et al. (2000Heger et al. ( , 2005 is often criticised on theoretical grounds as the various processes are not independent of one another (e.g., recently Chang & Garaud 2021, and references therein). For example, the dynamical and secular shears act on different timescales, and therefore their mutual use is physically contradictory. Maeder et al. (2013) proposed a diffusion coefficient accounting for the interactions between the different physical processes. While several studies have scrutinised these instabilities and resulting diffusion coefficients (e.g., Caleo et al. 2016;Goldstein et al. 2019;Barker et al. 2019Barker et al. , 2020Chang & Garaud 2021;Park et al. 2021, and references therein), a unified description of instabilities in rotating stars is still not fully complete.
In fully diffusive approaches as described above, two arbitrary scaling factors c and , introduced by Pinsonneault et al. (1989), are commonly adopted. If chemical gradients ∇ (Equation 14) develop, they may inhibit the efficiency of mixing. This is primarily due to ∇ serving as a stability criterion for the development of rotational instabilities (see e.g., Maeder 1997) since it appears directly in several of the individual diffusion coefficients used in Equation 16. To alter the effect of chemical gradients on mixing, the scaling factor is introduced such that ∇ is replaced by · ∇ when calculating stability criteria for various instabilities.
The parameter c , multiplying all individual diffusion coefficients in Equation 16, was first calibrated to c = 0.046 by Pinsonneault et al. (1989). This reduction in the efficiency of chemical mixing (compared to angular momentum transport) was needed to explain the observed lithium depletion in the Sun. However, recent studies (e.g., Prat et al. 2016) found that, at least, for the shear instability, both chemical mixing and angular momentum transport should have similar efficiencies when using the same diffusion coefficient. Heger et al. (2000) found that c = 0.033 with = 0.05 (which were the default options until recently) best reproduce the observed nitrogen enrichment in the 10 − 20 M mass range at Solar metallicity 17 . Yoon et al. (2006) concluded that when the angular momentum transport is very efficient (by using the magnetic term ST accounting for the Tayler-Spruit dynamo), then c = 0.033 should be used with = 0.1 instead of = 0.05. Using similar physical assumptions as Heger et al. (2000) and Yoon et al. (2006), Brott et al. (2011a) calibrated c = 0.0228 for a 13 M model based on the surface enrichment of early B stars in the LMC 18 and adopted = 0.1 from Yoon et al. (2006). Recently, Markova et al. (2018) found that this calibration produces insufficient mixing for more massive stars to be compatible with observations. Some subsequent modelling approaches even adopt a mixing efficiency parameter c that is a factor of 10 higher (Aguilera-Dena et al. 2020).
When using similar physics (assuming a purely diffusive equation to model angular momentum transport), Chieffi & Limongi (2013) obtained calibrations for c = 0.07 with = 0.03 and c = 0.2 with = 1.0 (correctly noting the degeneracy between these parameters). They also performed calibrations with different physics (using the advecto-diffusive equation of angular momentum transport) which yielded c = 1 with = 0.03 for chemical mixing. Recently, also using a physical approach different from the above mentioned ones, Costa et al. (2019) used intermediate-mass binary systems and constrained c = 0.17 with = 0.47. To be able to compare to previous works which used the same physics, we adopt c = 0.033 and = 0.1 (am_D_mix_factor = 0.033, am_gradmu_factor = 0.1) when using the Mix1 scheme 19 .

Mix2 scheme
Another commonly used mixing scheme ("Mix2" hereafter) was developed by Zahn (1992), Chaboyer & Zahn (1992), Maeder (1997); Maeder & Zahn (1998), Maeder & Meynet (2000). This scheme has been applied in the Geneva stellar evolution code ( , Eggenberger et al. 2008;Ekström et al. 2012;Georgy et al. 2013;Meynet et al. 2013;Groh et al. 2019;Keszthelyi et al. 2019;Murphy et al. 2021), as well as in modelling approaches using the (Potter et al. 2012a,b) and codes (Chieffi & Limongi 2013). Here, we adopt it in , which treats angular momentum transport in a fully diffusive scheme, unlike the above mentioned approaches. Therefore, a direct comparison to previous works is not possible. 17 Heger et al. 2000 also comment that for an initial rotation of 200 km s −1 and fixed = 0.033, ≥ 0.25 is inconsistent with observations in the 30-60 M mass range. 18 These values were also adopted for their Solar and SMC models. 19 We note that presently there is a growing amount of evidence that such a reduction in the efficiency of chemical mixing caused by hydrodynamical instabilities is likely not needed at all. Instead, there exist other processes that are simply more efficient in transporting angular momentum than chemical elements, the prime candidates being internal gravity waves and internal magnetic fields (e.g., Aerts et al. 2019, and references therein).
The major difference is that given the efficient diffusive angular momentum transport, strong shear mixing cannot develop. Consequently, in our models with the Mix2 scheme, the main chemical element transport is via meridional currents during most of the main sequence evolution. This is not the case in the models of Ekström et al. (2012), where the advective treatment of angular momentum transport allows for shears, which may also become the dominant process of transporting chemical elements. As we will see (Section 3.1, Section E), the Mix2 scheme leads to quasi-chemically homogeneous evolution for the entire main sequence of our nonmagnetic models. Since such a behaviour is expected to be rare, we may consider the adaptation of this mixing scheme in our models as a limiting case for very efficient mixing. The effective diffusion coefficient for chemical mixing combines the effects of meridional currents and horizontal turbulence, where the radial component of the meridional circulation is with the pressure, the density, the gravitational acceleration, the temperature, the luminosity, and terms which depend on the distribution of angular velocity and mean molecular weight 20 , and the ratio of the variation of the density to the mean density m . The horizontal turbulence is adopted as: where h is a constant set to unity (see Chaboyer & Zahn 1992), and ( ) expresses the radial dependence of the horizontal component of the meridional circulation. The horizontal component is expressed as ( ) 2(cos Θ), where 2 is the second Legendre polynomial and Θ is the co-latitude. We set ( ) = ( ) as a reasonable approximation. Then, The diffusion coefficient accounting for vertical shear mixing is derived by Maeder (1997) as: where energ is a free parameter set to unity. is the local pressure scale height, the gravitational acceleration, Ω the angular velocity, and the radius. Finally, in the Mix2 scheme, the diffusion coefficient applied in Equation 12 is: 20 The full expression of these terms is given by Maeder & Zahn (1998). In our approach, we simplify this expression and adopt only the leading term which is the first term of ★ Ω as described by Maeder & Zahn (1998). Since it is a smaller term, we set to zero.
In this case, the free parameters c and are not used in calculations. Consequently, we do not apply them in our Mix2 model calculations either (am_D_mix_factor = 1, am_gradmu_factor = 1). To our knowledge this mixing scheme is implemented in for the first time. The sum of diffusion coefficients, dominated by meridional currents ( eff , Equation 17) in the solid-body rotating case are comparable in shape to the default approach which uses ES as derived by Kippenhahn (1974) and Pinsonneault et al. (1989). However, the amplitudes are not equal (as shown in Figure 3).
We note here that there is some confusion in the literature regarding the work of Chaboyer & Zahn (1992). Heger et al. (2000) (and following publications) state that Chaboyer & Zahn (1992) found c = 1/30 based on a theoretical approach. The work of Chaboyer & Zahn (1992) does not introduce any scaling factors. eff , describing the transport resulting from the interaction between meridional currents and the strong horizontal turbulence is obtained by integrating the equation for the transport of chemical elements over latitudes. This integration gives rise to the numerical term of 1/30 (in their Equation 16 and our Equation 17), resulting from the decomposition of the meridional velocity in Legendre polynomials. This is not a scaling factor to match observations and it does not apply to any other diffusion coefficient. Similarly, in Equation 9 the numerical term 2/3 is not an arbitrary scaling factor that one would tailor to observations.

RESULTS
In this section, we first consider non-magnetic models, and then a fiducial model with ini = 20 M , Ω ini /Ω crit,ini = 0.5, ini = 0.014, and eq,ini = 3 kG within the INT/Mix1 scheme, and follow changes in its stellar structure and evolution. In particular, we will first vary the mixing and braking schemes to investigate the impact on stellar structure models in Section 3.1.1 and abundances in Section 3.1.2 and in evolutionary models in Section 3.2.1. Then, a typical HRD evolution of the INT and SURF magnetic models in the full mass range (3 -60 M ) will be addressed in Section 3.2.2, followed by predictions for the Kiel diagram in Section 3.2.3. Nitrogen abundances and other schemes are also shown in Appendix E. Finally, the initial magnetic field strength and metallicity will be varied within the evolutionary models in Sections 3.2.5 and 3.2.6. Figures 2, 3, and 4 show 20 M models at solar metallicity. These structure models are at half-way through their core hydrogen burning phase 21 (defined as core / core,init = 0.5). In Figure 2, we show non-magnetic models in the Mix1 (top panel) and Mix2 (lower panel) chemical mixing schemes. (Since magnetic braking is not applied, there is no braking scheme in these cases, hence we refer to these models as "NOMAG".) In Figure 3 Figure 4. Note that the models with the adopted braking and mixing schemes correspond to different stellar ages when the core 21 The ZAMS and TAMS structure models of the four schemes are shown in Figures D1 -D4 in the Appendix. hydrogen is half-way depleted in each model (indicated in the title of the panels), given the different evolution resulting from the change in physical assumptions.

Chemical mixing and angular momentum transport
When fossil magnetic fields are not considered in the models (NOMAG case; eq = 0), the Mix1 and Mix2 chemical mixing schemes produce drastically different results. In the Mix1 scheme, the mean molecular weight (shown with magenta line and on the right ordinate) drops rapidly at the convective core boundary (the zones with convective overshooting are shown with orange in Figure 2). Near this region the GSF instability dominates (see, e.g., the recent studies of Caleo et al. 2016;Barker et al. 2019Barker et al. , 2020Chang & Garaud 2021). Once the chemical composition is stabilised, meridional circulation drives chemical mixing (red line) and angular momentum transport (dotted line). Note that the assumption of c = 0.033 reduces the efficiency of all instabilities considered for chemical mixing compared to the efficiency of the same instabilities used for angular momentum transport. The angular velocity profile (right panel) remains completely flat in the stellar envelope, with a small break at the core boundary. Thus the model is very close to solid-body rotation. The NOMAG/Mix2 model reveals a very efficient mixing, with an almost flat mean molecular weight profile indicating close-to chemically homogeneous evolution. The dominant transport is via the effective diffusion coefficient (Equation 17). Importantly, the diffusion coefficients at the core-envelope boundary are smooth. This is a critical region that allows for mixing up material from the core to the surface. Note the much larger convective core (in grey) compared to the Mix1 model. The specific angular momentum (blue line, right panel) and angular velocity profiles are smooth throughout the star. Ω slightly decreases near the surface as a result of mass loss, however, this model is also very close to solid-body rotation.
In the INT braking scheme ( Figure 3) the angular velocity profile is completely flat. The star is rigidly rotating due to the assumed high diffusivity for transporting angular momentum attributed to the magnetic field, albeit the angular rotation is much lower than in the NOMAG model due to magnetic braking (uniformly) lowering the specific angular momentum. The rigid rotation does not allow shears to develop and transport angular momentum or chemical elements. Therefore in these models the chemical enrichment is entirely driven by meridional currents. In the "standard" description (Mix1 scheme) a gap in the transport develops above the overshooting region, corresponding to steep chemical gradients, as seen from the large drop of the mean molecular weight (magenta line, right ordinate) at the core boundary. Despite the mitigating effect of = 0.1 in these models, the inefficient mixing above the core boundary will prevent a very efficient surface enrichment and overall mixing inside the star. If the gap existed throughout the entire early evolution, it would completely inhibit surface enrichment. However, the gap is not present initially -see top panels of Figures D1 and D2 -, when the mixing and corresponding enrichment are prominent. With internal magnetic braking, all layers of the star lose angular momentum, therefore the shape of the specific angular momentum profile remains unchanged whereas its overall value decreases over time. In the Mix2 scheme, eff never becomes zero close to the convective core. This allows for a smoother composition gradient and more overall mixing, therefore differences in surface abundances are expected. On the other hand, the specific angular momenta are not so different between the Mix1 and Mix2 schemes in the INT models. Ω is smaller in the Mix2 model, but this quantity also depends on the radius of the star. Since the model in the INT/Mix2 scheme takes more time than the INT/Mix1 to deplete hydrogen in its core due to the more efficient chemical mixing, at half-way through its core burning stage it has a larger radius. We also note that the NOMAG/Mix2 model produces a more efficient mixing than the INT/Mix2 model. As a consequence, the INT/Mix2 model results in a smaller convective core than the non-magnetic case.
In the SURF braking scheme (Figure 4), one major difference is that there is no overall solid-body rotation. Let us recall that in the SURF models only the outer layers enclosing the top 20% of the total mass of the star are assumed to lose angular momentum (see Equation 10) and have an increased diffusivity for angular momentum transport (via Equation 11). In a 1D diffusive scheme, angular momentum flows from inner to outer layers. The angular velocity profile is flat in those layers where the diffusivity is increased. Depending on the mixing scheme, the composition and also the mixing processes are rather different. The magnitude and radial dependence of the diffusion coefficient for angular momentum transport also determines the angular velocity in the rest of the stellar envelope. We see that in the SURF/Mix1 model the angular velocity is also roughly constant between 0.6 and 0.8 / ★ , and it gradually changes closer to the boundary of the stellar core, as it is the case for AM (left panel, dotted line). In the SURF/Mix2 model, AM changes more abruptly at around 0.8 / ★ , whereas it is roughly constant again closer to the stellar core. Braking for a given magnetic field strength is less efficient in the SURF scheme than in the INT scheme since the Alfvén radius is smaller for a quadrupole field than for a dipole field as defined by Equations 3-5. Given the shape of the specific angular momentum profile (right panel, blue line), the SUFR/Mix1 model has a break closer to the stellar core, while the SURF/Mix2 model has a break closer to 0.8 / ★ . This means that the SURF/Mix2 model can more easily exhaust its surface reservoir of angular momentum.
Certainly, further research is required to investigate how angular momentum transport and magnetic braking work for more complex magnetic field configurations and how they could be implemented in 1D stellar evolution models. Overall, the results from the SURF approach may be considered similar to the works of Meynet et al. (2011) and Paper I, where magnetic braking was only applied to the uppermost stellar layer.
In the SURF/Mix1 scheme, the GSF instability can efficiently transport chemical elements near the core boundary. This instability acts on a dynamical timescale and therefore can vary from timestep to timestep. In the upper envelope meridional currents remain efficient. In the Mix2 scheme, the free parameters controlling mixing efficiency are not applied ( = 1, c = 1), and the SURF/Mix2 scheme is thus the most efficient in chemical mixing. This is also evidenced by the larger convective core size compared to the three models in the other magnetic schemes. In fact, the convective core size of the SURF/Mix2 model is similar to that of the NOMAG/Mix2 model. Strong gradients of chemical elements do  not develop near the core boundary. Shear mixing remains efficient in the entire envelope to transport chemical elements. Figure 5 shows the abundances of helium, carbon, nitrogen, and oxygen. For three models (INT/Mix1, INT/Mix2, SURF/Mix1) the convective hydrogen core sizes are comparable. The SURF/Mix2 scheme, which is the most efficient in mixing, leads to a much larger core. In this model, the average helium content in the stellar envelope is much higher and the surface will also show this increased abundance already on the main sequence. Carbon is slightly depleted during the CNO-cycle (it becomes most depleted in a thin layer close to the core boundary), however the surface carbon abundance is minimally changed in the first three models. In contrast, the SURF/Mix2 model produces an almost uniform distribution of carbon inside the star, it is completely mixed in the envelope without any gradients. Nitrogen excess is produced during the CNO-cycle, and therefore, its surface abundance is a crucial measurement to infer the efficiency of internal mixing. All models produce a surface nitrogen enrichment, except the INT/Mix1 scheme. Here a strong gradient develops between the core and the surface. While a coresurface gradient is also present in the INT/Mix2 and SURF/Mix1 models, their envelopes have a somewhat higher mean nitrogen abundance and their surfaces are slightly enriched in nitrogen. The SURF/Mix2 model has an almost homogeneous nitrogen distribution in its envelope. Oxygen is depleted during the CNO-cylce. Similar to carbon, oxygen can still remain abundant in the envelope in the first three models. The SURF/Mix2 model yields a close-tohomogeneous oxygen distribution throughout the star. In contrast to carbon, in all four models oxygen is depleted in the entire stellar core. For example, the core to surface oxygen abundance can differ by an order of magnitude in the INT/Mix1 model.

Abundances of He, C, N, O
The reason why strong gradients can develop (and remain in most models) near the core boundary is related to the drop in chemical mixing, identifiable by drops and gaps in the diffusion coefficients (c.f. Figures 3-4), which in turn depend on the composition gradients. The velocity of the diffusion is zero when there is no composition gradient and it increases when the gradient increases. This is another effect that leads to reducing the diffusion coefficient. For helium, the difference between the core and the envelope grows gently, while for nitrogen it grows faster because nitrogen is enhanced very rapidly in the core. Thus for a given diffusion coefficient, nitrogen will diffuse more rapidly than helium. In the Mix2 scheme, the key differences between the INT and SURF models result from their different angular velocity profiles. The INT models lose angular momentum in all layers and thus mixing becomes less efficient overall. The SURF/Mix2 model has a more massive convective core at the same evolutionary stage (c.f. Figure 4), and the envelope closely reflects on the core composition as this model has an almost homogeneous distribution of chemical elements. In both INT/Mix1 and INT/Mix2 schemes, envelope mixing has a similar overall efficiency. However, in the INT/Mix2 scheme, the near core mixing is more efficient (see above), greatly impacting the measurable surface abundances of carbon, nitrogen and oxygen. The surface abundances, especially of nitrogen, are sensitive to the chosen braking and mixing schemes, especially with the Mix2 scheme reflecting more closely the core composition than the Mix1 scheme.

Impact of magnetic braking and chemical mixing schemes
for 20 M models Figure 6 shows the fiducial 20 M model (INT/Mix1) as well as initially identical models but in the other 3 schemes in the HRD, Kiel diagram, and Hunter diagram 22 (Hunter et al. 2008(Hunter et al. , 2009). Here we address the impact of using the different braking and mixing schemes for otherwise identical models with the same initial mass and initial magnetic field strength.
The models within the Mix1 chemical mixing scheme result in closely overlapping tracks on the HRD and Kiel diagrams. However, as evidenced from the Hunter diagram (right panel of Figure 6), the spin-down and chemical enrichment are different when considering the different magnetic braking schemes. The 20 M INT/Mix1 model (solid line) produces essentially no observable surface nitrogen enrichment in this configuration. The SURF/Mix1 model (dashed-dotted line), on the other hand, maintains a higher angular velocity in the inner regions of the star by braking only the upper layers (see Section 3.1). This is why mixing remains more efficient and a larger amount of nitrogen is mixed to the stellar surface (0.45 dex). Model predictions within the Mix2 chemical mixing scheme produce a much more efficient chemical mixing than the Mix1 scheme. However, the differences between the INT/Mix2 and INT/Mix1 schemes are relatively modest on the HRD and Kiel diagrams, while the Hunter diagram shows large deviations. The INT/Mix2 model reaches a surface nitrogen abundance that is about 0.4 dex higher than the baseline value.
Contrary to the first three cases, the SURF/Mix2 model has an extended blueward evolution, shown on the HRD and Kiel diagrams. Such a feature is commonly associated with blue stragglers, merger products, and quasi-chemically homogeneous evolution (e.g., Maeder 1987;Yoon et al. 2006;Knigge et al. 2009). Given the very efficient mixing in this model, it is nearly chemically homogeneous (c.f. Figure 5). Indeed, the SURF/Mix2 model produces the highest nitrogen enrichment, more than 1 dex compared to the baseline in this configuration.
Even without considering magnetic braking, the two mixing schemes are considerably different, which leads to different evolutionary tracks. The assumptions regarding the mutual effect of the magnetic field and chemical mixing allow for a range of behaviours (see also Paper III). Models within the Mix1 scheme have modest differences as a result of using the different braking schemes (INT/SURF). Models within the Mix2 scheme produce very efficient mixing, which is more easily quenched in the INT models than in the SURF models. The INT braking scheme, in contrast to the SURF, decreases the overall rotation rate, and the INT/Mix2 models do not evolve significantly blueward on the HRD. However, quasichemically homogeneous evolution is achieved in the case of nonmagnetic Mix2 models, as well as in some magnetic SURF/Mix2 models depending on the initial field strength and stellar mass (see Figs. E3 and E4, in the Appendix). For example, the magnetic model in the SURF/Mix2 scheme leads to quasi-chemically homogeneous evolution for most of the main sequence (≈ 9 Myr out of 12 Myr) of a 20 M model at solar metallicity with a 3 kG magnetic field. How- 22 We use the spectroscopic definition of nitrogen abundance for the Hunter diagram, which is log( / ) + 12, where and are the surface number fraction of nitrogen and hydrogen. Note that mass fractions (which are the typical output quantities from evolutionary grids) need to be translated to number fraction by appropriate scaling.  ever, for the same initial field strength, the initially 60 M model only experiences a brief phase of quasi-chemically homogeneous evolution and then evolves redwards ( Figure E4). In summary, the braking and mixing schemes can drastically change the main observable characteristics. For example, the 8 Myr isochrone spans over a 10 kK effective temperature range (left panel of Figure 6) despite the models being initially completely identical. (When the "extreme" case of SURF/Mix2, which is assumed to represent stars with complex magnetic fields and very efficient mixing, is not considered the difference is less, around 2 kK.) Thus the uncertainties associated with braking and mixing schemes in evolutionary model predictions are significant. From the Hunter diagram we can conclude that various braking and mixing schemes can cover a wide range of rotation rates and surface nitrogen abundances. Three models already reach slow rotation (< 50 km s −1 ) with high surface gravities. Both INT models take less than 6 Myr to achieve this, while it is slightly over 8 Myr for the SURF/Mix2 model. The SURF/Mix1 model reaches the TAMS with a somewhat higher rotation rate than the other models. E2 and E4 in the Appendix for the SURF schemes) and we demonstrate the impact of a magnetic field with an initial equatorial field strength of 500 G.

HRD evolution of a grid of magnetic models
In the INT models, strong magnetic braking leads to a rapid decrease of surface rotational velocity in the entire mass range from 3 to 60 M . This implies that rapidly-rotating (single) magnetic massive stars are expected to be young and close to the ZAMS on the HRD. Some quantitative differences arise from the assumptions of the chemical mixing schemes. Nonetheless, within the INT magnetic braking scheme these differences are small on the HRD, albeit it could affect the parameter determination (current age and mass) of known magnetic stars. Figure 8 shows the Kiel diagram. In the SURF models, radial differential rotation develops between the stellar core 23 and surface, as indicated with the grey contour lines. This is because we assume that a magnetic field with a complex geometry would only exert a strong torque on the near-surface layers by considering that organised, strong magnetic flux is not present in deeper stellar layers (see also Braithwaite & Nordlund 2006;Braithwaite 2008). Thus angular momentum transport in the deepest layers of the stellar envelope is dominated by -less efficient -hydrodynamical instabilities. On the other hand, in the INT models complete solid-body rotation is maintained throughout the main sequence. Similar to Figure 6, where the 20 M models are discussed, we see that the Mix1 and Mix2 chemical mixing schemes (within the SURF braking scheme here) result in notable differences in observable stellar parameters, which are the most prominent in the 5-10 M range.

Differential rotation on the Kiel diagram
If the assumed magnetic field geometry remained unchanged throughout the main-sequence evolution, the prediction for magnetic massive stars is that differential rotation could be best identified in hotter, more massive, and more evolved (lower log ) stars. An important caveat nevertheless is the time evolution of complex magnetic fields. According to Braithwaite (2008), complex fields 23 Here, we use the output quantity of core angular velocity "center_omega" to determine the rotation of the stellar core. The exact radius to obtain Ω c is not essential since the entire stellar core has the same angular velocity. may simplify to a dipolar form since higher-order harmonics diffuse more rapidly. Indeed, Shultz et al. (2019b) finds that evolved magnetic stars tend to have simpler geometries (however, see also Kochukhov et al. 2019). Insofar it remains unclear what the exact diffusion time of complex magnetic flux tubes would be (perhaps still longer than the main-sequence lifetime) and whether the corresponding effects which are modelled in evolutionary codes (angular momentum transport, magnetic braking, mass-loss quenching) would change significantly. We evaluated the models with different initial field strengths and found that a stronger magnetic field is able to achieve a higher degree of differential rotation in comparison to models with lower field strengths. It would therefore be of great importance to obtain seismic data of a sample of magnetic massive stars.

Mass-dependent rotational evolution
Both the HRD and Kiel diagrams presented above 24 show the same distinctive feature. Namely, irrespective of how fast the spin-down per given scheme is, the mid-mass range models (≈ 5-10 M , the transition typically taking place at around 5-7 M depending on model assumptions) always maintain a higher rotational velocity than models with other masses. This is the most striking on the left panel of Figure 8 for the SURF/Mix1 scheme. (The exact model behaviour was also recognised for a 5 and 10 M model in Paper II.) The distinctive feature is a consequence of the spin-down of the models, which does not scale linearly with mass. For higher-mass models, wind mass loss contributes to the spin-down. The spindown time to reach a given surface rotation takes a lower fractional age when the mass increases for initial masses higher than 10 M . Since stars below about 5 M have much longer nuclear timescales than higher-mass stars, even if magnetic braking is less efficient Figure 8. Kiel diagram of magnetic evolutionary models with eq,ini = 3 kG, Ω ini /Ω ini,crit = 0.5 at solar metallicity ini = 0.014. within the SURF/Mix1 (left) and SURF/Mix2 (right) schemes. Grey contour lines indicate the unitless degree of differential rotation, quantified as the ratio of core to surface angular velocity. Note that Figure 7 shows INT models where solid-body rotation (Ω core /Ω surf = 1) is achieved throughout the main sequence (see also Figure D1). The initial masses of the models decrease from left to right. The colour-coding shows the surface equatorial rotational velocity. due to weaker winds, the available timescale allows for braking the rotation at a lower fractional main-sequence age. This results in stars below 5 M rotating more slowly at the TAMS than in the 5-10 M range. We emphasise that the extent of the tracks on the HRD is not linear with stellar age: a significant time evolution can take place in a narrow location on the HRD close to the ZAMS as demonstrated by the isochrones. For example, the 50 M model spends the first 2 Myr of its evolution while decreasing its eff by only about 5 kK, whereas in the second 2 Myr of its evolution, its eff decreases by about 15 kK ( Figure 7). Moreover, from 0 to 4 Myr, the 25 M track evolves roughly 0.2 dex in luminosity and a few kK in eff , which is a typical range of observational uncertainties depending on data quality and knowledge of distance and extinction. The precise age determination of young stars especially at initial masses below 25 M becomes rather challenging with uncertainties well exceeding 1 Myr.
In summary, stars evolve slowly in the HRD at the beginning of the main sequence phase, and thus suffer strong magnetic braking while not evolving away from the ZAMS (in effective temperature and luminosity). However the interplay between the dependence on the initial mass of the meridional current velocity, the evolutionary timescale near the ZAMS and the evolution of the radius produces a small bump of the surface rotational velocity in the initial mass range of 5-10 M (see also Paper II). At the TAMS, this gives rise to a slower rotational velocity in models below 5 M than models in the mid-mass range. Figure 9 shows the impact of the initial equatorial magnetic field strength within the INT/Mix1 scheme on the HRD, Kiel, and Hunter diagrams for models from 3 to 60 M . The colour-coding of the surface rotational velocity on the HRD and Kiel diagrams show that for initial masses above 30 M , stellar winds play a significant role in depleting the angular momentum reservoir, and thus even without magnetic fields those stars can significantly spin down on the main sequence. However, in the mass range from 3 to 30 M single-star models would not undergo a dramatic angular momentum loss unless they were strongly magnetised. As demonstrated in the figure, the stronger the magnetic field, the more rapidly the surface rotation brakes. In particular, already a 3 kG equatorial field would produce a (sub)population of stars whose surface rotation is less than 50 km s −1 throughout essentially the entire main sequence.

Impact of varying the initial equatorial magnetic field strength
The Kiel diagram shows yet another consequence of magnetic fields. Apart from the highest-mass models (> 30M ), the models with stronger magnetic fields tend to reach the end of the main sequence with higher surface gravities. For example, for a 10 M model, the TAMS value of log increases from 3.3 to 3.4 and to 3.5 from the 0 kG to the 0.5 kG and to the 3 kG models, respectively. Non-magnetic models maintain a higher rotation and thus the mixing remains more efficient. This mixing (if not too strong to keep the star in a bluer position in the HRD) tends to enlarge the convective core and consequently, extend the width of the main sequence towards lower effective temperatures and higher surface gravities.
The Hunter diagram reveals that the magnetic models may strongly deviate from non-magnetic model predictions. The stronger the magnetic field, the more rapidly rotation brakes, and the less nitrogen can be mixed to the stellar surface. We identify and demonstrate that there exists a cutoff magnetic field strength, above which no surface enrichment is expected given that it leads to a shorter magnetic braking timescale compared to the rotational mixing timescale. However, this cutoff field strength is strongly model and parameter dependent, and thus the exact value varies for given stellar mass, initial rotation, metallicity, and mixing scheme, amongst others. We evaluate this in Section 4.2.  Figure 10 shows the impact of the initial metallicity (here, for LMC and SMC values; the Solar metallicity models with the same input are shown in the lower panels of Figure 9) for 3 to 60 M models with initial equatorial magnetic field strengths of 3 kG within the INT/Mix1 scheme on the HRD, Kiel, and Hunter diagrams. At lower metallicity, the ZAMS is shifted to higher effective temperatures given that the stellar models are more compact due to the lower opacity and lower mean molecular weight. Specifically, the lower CNO abundances impose further contraction of a star to initiate core burning.

Impact of metallicity
Magnetic braking, in our formalism, is metallicity independent. However, rotational mixing is not. The various mixing prescriptions depend on chemical composition and their gradients, which in turn affects the evolution of surface rotational velocity as most prominently revealed on the Kiel diagrams. For example, for a given value of a diffusion coefficient, the mixing timescale is ≈ 2 / chem . Since stars are more compact in lower , the timescale becomes shorter. Consequently, the changes in rotation and surface abundances can be more impacted in lower metallicity stars.
Similarly, the highest relative nitrogen enrichment is seen when metallicity is the lowest (see also e.g. Brott et al. 2011a;Georgy et al. 2013). Let us recall that Figure 10 shows the INT/Mix1 models, which -in our approach -are the lowest estimates for the surface enrichment (c.f. Figure 6). The other schemes predict higher surface nitrogen enrichment when combining the effects with low metallicity. The trends produced in magnetic massive star models are unique since they lead to simultaneous surface nitrogen enrichment and a rapid spin-down of the stellar surface (c.f. Paper I, Paper II). The rapid spin-down could otherwise only be expected for very massive stars with extremely strong stellar winds.

Slowly-rotating nitrogen enriched stars in the LMC
The projected rotational velocities of massive stars in the Magellanic Clouds appear to follow a bi-modal distribution (e.g., Dufton et al. 2013Dufton et al. , 2018Dufton et al. , 2019Dufton et al. , 2020Rivero González et al. 2012;Ramírez-Agudelo et al. 2013. The bi-modality is also observed for intermediate-mass stars up to about 5 M (e.g., Bastian et al. 2020;Sun et al. 2021). The observed slowly-rotating red main sequence stars and rapidly-rotating blue main sequence stars are thought to be evidence for main sequence splitting (e.g., Bastian et al. 2020) and an extended main sequence turn-off (e.g., D'Antona et al. 2015). It has been suggested that the low-velocity peak might be caused by magnetic braking (e.g., Wolff et al. 1982;Sun et al. 2021). Shultz et al. (2018) demonstrated that the dichotomy in sin between Galactic B-type stars with and without magnetic fields is at least qualitatively consistent with the lower sin values observed in the magnetic population. For observed massive stars in the Magellanic Clouds, a notable fraction of slow-rotators were found to show measurable nitrogen enrichment, which challenges typical, nonmagnetic single-star evolutionary models (Lennon et al. 1996(Lennon et al. , 2003Dufton et al. 2006Dufton et al. , 2013Dufton et al. , 2018Dufton et al. , 2019Dufton et al. , 2020Rivero González et al. 2012;Ramírez-Agudelo et al. 2013McEvoy et al. 2015;Grin et al. 2017).
The nitrogen-enriched slow-rotators (also known as "Group 2" stars, Hunter et al. 2008) correspond to roughly 20% of the population in the LMC (e.g., Hunter et al. 2008;Brott et al. 2011b;Grin et al. 2017;Dufton et al. 2018). In the Galaxy, the observed incidence rate of fossil magnetism is found to be ≈ 10% (e.g., Fossati et al. 2016;Grunhut et al. 2017;Sikora et al. 2019a) and it has previously been suggested that at least some of the Group 2 stars could be explained by magnetism (Meynet et al. 2011;Potter et al. 2012b, Paper I). This would require an incidence rate of fossil Figure 11. Hunter diagram of magnetic single-star evolutionary models with eq,ini = 3 kG, Ω ini /Ω ini,crit = 0.5 at LMC ( ini = 0.0064) metallicity within the four schemes. Models within the SURF/Mix1 scheme are also shown with white lines since they overlap with the INT/Mix2 models. Additionally, two sets of non-magnetic (NOMAG) models are shown within the Mix1 (grey) and Mix2 (black) schemes. For visualisation purposes, we reduced the numerical noise in the latter case. Models with initial masses from 15 to 60 M are shown. The actual surface equatorial rotational velocity of the models is scaled by sin( /4) to account for an average inclination angle. The coloured area corresponds to our definition of Group 2 stars. The colour-coding of the models shows the logarithmic surface gravity. Observations are shown with circles and squares, respectively. A typical reported uncertainty in the observed nitrogen abundances is about 0.1 dex. magnetism in the LMC that is likely higher than the 10% observed in the Galaxy 25 . In addition, the (initial) magnetic field strength distribution is not yet known in our galaxy or in other metallicity environments; however, see Petit et al. (2019) and Cerrahoğlu et al. (2020) for theoretical models. It could also be that the Group 2 stars require an additional channel to explain all observations. In fact, binarity has been suggested by, e.g., Song et al. (2018b). Figure 11 shows the Hunter diagram with models representative of LMC metallicity ( = 0.0064). Let us recall that we 25 For example, lower metallicity environments might favour a higher incidence rate of stars with fossil fields if the convective expulsion scenario, due to the sub-surface iron opacity bump, regulates the incidence rate (Jermyn & Cantiello 2020). specifically adopted an initial nitrogen abundance in our models of log( / ) + 12 = 7.15 from Dopita et al. (2019) to produce evolutionary models guided by available empirical baseline abundance determinations. Here we demonstrate some of the complex parameter-space dependences that magnetic single-star models produce, albeit strongly depending on the model assumptions, especially the mixing and braking schemes (see also Meynet et al. 2011;Potter et al. 2012b, Paper I, Paper III). We display the non-magnetic models as well as magnetic models in the two braking and two mixing schemes. The initial equatorial magnetic field strength is 3 kG and the initial rotation is set by Ω ini /Ω ini,crit = 0.5. These assumptions produce evolutionary models, which over time reasonably approximate mean values measured from observations (magnetic field strengths from, e.g., Shultz et al. 2018 for Galactic magnetic B-type stars, and rotational velocities from, e.g., Dufton et al. 2013 for massive stars in the Magellanic Clouds). For example, the 3 kG initial ZAMS magnetic field strength weakens by roughly an order of magnitude (since eq ∝ −2 ★ over time) at the TAMS to 300 G. Only models with initial masses from 15 to 60 M are shown given that the available instrumentation allows magnitude limited observation of bright LMC stars that are more massive than ≈ 15 M (e.g., Schneider et al. 2018). Thus the models are the most representative of O (and early B-type) stars.
Non-magnetic models (black and grey lines for the Mix1 and Mix2 schemes, respectively) mostly show high rotational velocities. Close to the TAMS, the NOMAG/Mix2 models spin down efficiently and yield a high N/H ratio. However, this prediction is associated with producing Helium stars since these models experience Wolf-Rayet type mass loss due to their quasi-chemically homogeneous main sequence evolution (see Section D2). This helps to decrease the surface hydrogen abundance significantly, and hence the N/H ratio can further increase. Given the high effective temperatures that these models show close to the TMAS, we do not expect that the non-magnetic models should match the observations of typical Group 2 stars.
Chemical mixing remains challenging to constrain. The schemes that we assume in this work cover a large area on the Hunter diagram, which represents modelling uncertainties. In our models, for a given initial rotation and initial magnetic field strength, the INT/Mix1 and SURF/Mix2 models lead to the smallest and highest amount of nitrogen enrichment, respectively (see also Figure 6). The INT/Mix2 and SURF/Mix1 models produce similar results, although the latter models cover a narrower domain. In this sense, the Mix1 and Mix2 schemes are limiting cases in terms of the produced nitrogen enrichment from our models. In nature, the situation might be much more complex since, for example, the study of slowlypulsating B-type stars reveals a diverse range of mixing profiles (Pedersen et al. 2021). The variation of these profiles over time is not yet quantified.
The magnetic properties of stars in the LMC remain unknown. For this reason, the INT and SURF braking schemes represent assumptions and uncertainties that could only be resolved if, at least, upper limits on the magnetic field strengths were constrained. From the models we see that for a given initial magnetic field strength, the INT models produce less enrichment than the SURF models in a given mixing scheme.
In Figure 11 we show abundance measurements 26 of observed massive stars at the LMC made by Grin et al. (2017) and Dufton et al. (2018). Since the observations only allow derivation of the projected rotational velocity sin , we scale the actual rotational velocity in our models with sin( /4) to account for an average inclination angle. We only consider here observations with sin < 300 km s −1 .
The observations reveal a large scatter, which likely represents a range of initial conditions (mass, rotation rates, magnetic field strength, binarity, etc.) and current age. The data from Grin et al. (2017) and Dufton et al. (2018) indicate stars in their main sequence evolutionary stages with surface gravities systematically decreasing towards lower rotation rates. Our models show that magnetic braking typically yields slowly-rotating stars early on in the evolution, still with high surface gravities. Once the rotation is slow (and log is still ≈ 4.0), chemical mixing becomes inefficient and no further surface enrichment may be expected on the main sequence. This is the primary reason why magnetic models produce less surface enrichment than non-magnetic models with the same mixing assumptions (see also Figure E1). However, at this point, the magnetic models still evolve further in time, albeit their location does not change in the Hunter diagram. Thus the magnetic models decrease their surface gravities (to log ≈ 3.0) in a narrow region in the Hunter diagram, when their (assumed projected) surface rotational velocities are below 50 km s −1 in all cases, except for the SURF/Mix1 models where the spin-down is the least efficient (c.f. Section 3.2.1).
There are several other caveats, which hamper a quantitative comparison between the models and observations. For example, the mass determinations are uncertain and often rely on (non-magnetic) evolutionary models. As we demonstrate in Figure 11, even for an idealised situation where the braking and mixing schemes were known for a given magnetic field strength, the produced nitrogen enrichment is still a function of initial mass (see also Aerts et al. 2014;Maeder et al. 2014). In the INT/Mix1 scheme, the final nitrogen abundance becomes a factor of 2.5 higher when increasing the stellar mass from 15 to 60 M . This difference is less in the INT/Mix2 and SURF/Mix1 schemes (factor of ≈1.5), while this trend reverses for the SURF/Mix2 scheme.
Despite all these uncertainties, the models incorporating the effects of surface fossil magnetic fields can cover the region on the Hunter diagram where the "anomalous" Group 2 stars (slow rotation along with surface nitrogen enrichment, see Hunter et al. 2008) are located, which is not possible with standard main sequence evolutionary models of single stars (see Martins et al. 2017 for non-magnetic ≈ 30 M models in the Galaxy). In particular, for the slowly-rotating non-magnetic Mix2 models the produced nitrogen enrichment seems to be larger than indicated by observations, whereas the non-magnetic Mix1 models do not spin down sufficiently.
Since the parameter space is degenerate, not only the mixing and braking scheme could produce results that cover Group 2 stars but also the variation of initial magnetic field strength in a given scheme. A stronger initial field would yield less enrichment (Figure 9), possibly explaining the less nitrogen enriched stars, whereas a weaker initial field could be compatible with the most highly enriched stars (see Figure E1). To quantify this, we introduce and discuss the cutoff magnetic field strength in the next section.

Cutoff magnetic field strengths in the LMC
It is of interest to evaluate a critical value of the magnetic field that strongly impacts observable properties. In particular, what range of initial magnetic field strengths allow for producing Group 2 stars? To this extent, we define a cutoff (maximum) field strength max,N as the initial equatorial magnetic field strength in a given model that allows for producing more than 0.1 dex of surface nitrogen enrichment (in spectroscopic units) during its main sequence evolution. If the initial magnetic field strength is higher than max,N , then the mixing is inefficient due to the magnetic spin-down, and no nitrogen enrichment can be observed. Similar to the above definition, we may also define a cutoff (minimum) magnetic field strength that is the initial equatorial magnetic field strength needed to produce sufficiently slow rotators ( rot < 71 km s −1 ) by the end of the main sequence evolution. Any value higher than min,v will yield slowrotating models; however, values below min,v will still result in considerable rotational velocities at the TAMS. Figure 12 shows the range of possible initial equatorial magnetic field strengths that are Figure 12. Cutoff magnetic field strengths as a function of mass to produce surface nitrogen enrichment ( max,N , solid line) and slow rotation ( min,v , dashed line). Values above max,N will inhibit surface nitrogen enrichment, whereas values below min,v will not spin down the star sufficiently. Left/right panels show the INT/SURF models. The models are considered for initially Ω/Ω crit = 0.5 at LMC ( = 0.0064) metallicity. able to produce Group 2 stars given the constraints given above. max,N is shown with solid line and min,v with dashed line. The range is shown as a function of initial mass for the LMC (the Solar and SMC metallicity models are discussed in Appendix F). The cutoff field strengths depend on the initial rotation rates, metallicity, mass, and chemical mixing and magnetic braking schemes. We discuss now the latter three.
Models within the INT/Mix1 scheme result in the lowest value of max,N since this scheme is the least efficient in chemical mixing. However, stronger magnetic fields would brake the rotation faster than the timescale of rotational mixing. In contrast, models within the SURF/Mix2 scheme have such strong mixing that even a 50 kG equatorial field strength is insufficient to inhibit mixing 27 . A particular feature, a jump in min,v around 5 to 10 M , can be explained by the mass-dependent rotational behaviour of the models discussed in Section 3.2.4. In most models, the decreasing value of min,v for stars more massive than 10 M implies that stronger stellar winds aid the spin-down thus a weaker magnetic field is sufficient to achieve slow-rotating stars.
We find that equatorial magnetic fields of initially a few kG are able to produce Group 2 stars in the INT/Mix1 scheme. Nevertheless, the range of allowed initial field strengths is the most limited in this case. In fact, models below 6 M are in a "forbidden" range where the minimum field strengths needed to brake rotation are higher than the maximum field strengths allowed to produce nitrogen enrichment. In the INT/Mix2 scheme, a much larger range of initial field strengths are allowed to produce Group 2 stars, particularly from 22 M , where the lower limit drops to 250 G. For initial masses higher than 7 M , the upper limit is of the order of 10 kG to produce Group 2 stars.
Models within the SURF/Mix1 and SURF/Mix2 schemes (right panel of Figure 12) cover a wide range of possible initial field strengths. The SURF/Mix1 scheme produces a similar massdependent pattern as the models in the INT scheme. Namely, models from 10 M have systematically decreasing values of min,v . How-27 For visualisation purposes, we assigned 60 kG to those models where the maximum value in our grid of models (50 kG) was still insufficient to prevent nitrogen enrichment. ever, for initial masses lower than 17 M a 10 kG initial equatorial field strength is still needed to achieve slow rotation. Interestingly, for initial masses ≥ 19 M there is a dip in max,N , staying constant at 30 kG in contrast to the mass range of 7-18 M , where the 50 kG initial equatorial field strength still allows for producing nitrogen enrichment. The SURF/Mix2 models have a constant upper limit given by max,N , meaning that even the strongest initial magnetic field strength we considered in this study is not sufficient to prevent nitrogen enrichment on the main sequence. The lower limit given by min,v is constant for initial masses higher than 6 M . In general, we can conclude that the INT scheme favours lower values for the cutoff magnetic field strengths of max,N to produce Group 2 stars and the SURF scheme allows for higher values of max,N . It is quite remarkable that the upper limit remains roughly constant in the INT cases for stars more massive than about 10 M . In the SURF/Mix2 case, likely an unrealistically strong magnetic field would be needed to prevent nitrogen enrichment (as the 50 kG field is still insufficient). The Mix1 scheme tends to allow for a narrow range of values and Mix2 scheme covers a wide range of possible solutions.
In the LMC, both quantities used in our criteria, surface nitrogen abundance and (projected) rotational velocity, can be measured. Large-scale surveys dedicated to magnetic field measurements are not yet available in lack of high-resolution spectropolarimetry (however, see Bagnulo et al. 2017Bagnulo et al. , 2020. This means that our predictions can be used as constraints on the strengths of magnetic fields that might exist in slowly-rotating nitrogen-enriched ("Group 2") stars in the Magellanic Clouds. An initial equatorial magnetic field strength above min,v and below max,N will produce stars that can be identified as Group 2 stars.

Future work
The models presented in this work cover the main sequence phase of single stars. Logical extensions include calculating pre-main sequence models and continuing the computations to the post-main sequence phase to be able to scrutinise connections with end-products of stellar evolution, such as strongly magnetised white dwarfs and neutron stars (magnetars). The recently discovered link between magnetars and fast-radio bursts (CHIME/FRB Collaboration et al. 2020; Bochenek et al. 2020) further supports investigations of the magnetic field origin of magnetars (e.g., Spruit 2009;Makarenko et al. 2021). Our grid of models could be further extended to cover initial masses below 3 M and thus to compare with, for example, Ap stars. This requires some further considerations about the winds of these objects and the inclusion of atomic diffusion in the models. We assumed single-star models in this work. Some magnetic massive stars, such as Sco (Schneider et al. 2016(Schneider et al. , 2020Keszthelyi et al. 2021), may challenge this scenario. Nonetheless, the vast majority of OBA stars with fossil fields have characteristics that do not require invoking a merger event. In particular, fossil magnetic fields are detected in young massive stars, for example, in Ori E (Landstreet & Borra 1978;Townsend et al. 2010;Oksala et al. 2012;Song et al. 2021). Certain binary systems also present challenges to the merger scenario, such as the doubly-magnetic binary Lupi (Shultz et al. 2015) and the eclipsing late B-type binary HD 62658, which comprises two young nearly-identical stars in a circularised orbit, only one of which is magnetic (Shultz et al. 2019c). Thus single-star models presented in this work are a reasonable first approach; however, future work remains to address binarity and mergers in combination with fossil field effects. In particular, multiplicity is common among massive stars (Sana et al. 2012de Mink et al. 2013de Mink et al. , 2014 and while close magnetic binary systems are rare (e.g., Alecian et al. 2013;Shultz et al. 2018), the mutual impact of magnetism and tidal interactions need to be further studied (e.g., Song et al. 2018a;Vidal et al. 2018).
The models computed in this work can be confronted with observations of known magnetic massive stars. They will complement previous approaches which relied on grids of stellar evolution models that did not include surface magnetic field effects (e.g., Brott et al. 2011a;Ekström et al. 2012;Chieffi & Limongi 2013;Choi et al. 2016) to infer stellar parameters and ages of magnetic stars. The differences are expected to be most pronounced for higher-mass stars , whereas -within the framework considered here -the lower-mass, A-type stars should be less impacted (Deal et al. 2021). Nonetheless, the available TESS data and continuous spectropolarimetric monitoring can be used to constrain accurate rotation periods of such stars and directly compare with evolutionary models incorporating magnetic braking.
The INT models represent stars with strong, predominantly dipolar fields that are commonly identified in the sample of known magnetic massive stars. The SURF models are a limiting case motivated by stars with complex magnetic fields. In the future, the implementation of different magnetic field configurations and their time evolution could be considered to improve the present models.
The internal mixing efficiency remains uncertain in evolutionary modelling. In our models, quasi-chemically homogeneous evolution during the main sequence develops in the non-magnetic Mix2 case. In the INT/Mix2 models this behaviour is prevented by the efficient overall spin-down already when initially weak magnetic fields are considered, whereas in the SURF/Mix2 models a very efficient mixing still remains. The occurrence and duration of quasi-chemical evolution depends on the initial field strength and stellar mass; some SURF/Mix2 models evolve significantly bluewards, whereas some models turn to a redward evolution after mixing becomes less efficient (see Figure D4). Quasi-chemically homogeneous evolution is expected to be rare in nature; however, it may be a crucial channel, for example, for some supernova events, gamma-ray bursts, or gravitational wave sources (e.g., Georgy et al. 2012;Martins et al. 2013;de Mink & Mandel 2016;Szécsi 2017).
Observational studies should help constrain the mixing efficiency and ultimately the physical mixing processes. As such, it would be beneficial to further study our models and confront them with measurements of surface nitrogen abundances in magnetic massive stars (Morel et al. 2008(Morel et al. , 2015Aerts et al. 2014;Martins et al. 2012Martins et al. , 2015, as well as studies which identified anomalous trends on the Hunter diagram in the LMC and SMC (e.g., Dufton et al. 2020, and see Section 4.1).
Finally, our current understanding of magnetic field evolution is still incomplete and, in particular, how different field geometries evolve over time is largely unconstrained. It will therefore be valuable to explore various field evolution scenarios, for example, magnetic flux decay (Shultz et al. 2019b, Paper III). This might also lead to a time-dependent magnetic braking scheme depending on the relative dissipation timescales of various complex field components.

CONCLUSIONS
In this work we present the most extensive grid of stellar structure and evolution models taking into account the effects of surface fossil magnetic fields. The grid is publicly available on Zenodo and we recommend that, while acknowledging the uncertainties, it could be used to infer stellar parameters of known magnetic massive stars. No particular braking (INT/SURF) or mixing (Mix1/Mix2) scheme can be preferred at this time, although we do assume that if the field geometry is known, the INT case is applicable for dipolar fields and the SURF case for more complex geometries. Thus we consider the four schemes as limiting cases and, as we demonstrate in this work, the differences between these subgrids can substantially impact the determination of stellar parameters. It is therefore essential to confront the Mix1/Mix2 mixing schemes with spectroscopic studies even for stars where surface magnetic fields are not detected. In all cases (2 mixing schemes and 3 metallicities), we provide a subgrid of non-magnetic evolutionary models. Furthermore, the grid of models is suitable for population synthesis studies, which thus far have neglected magnetic field effects and magnetic massive stars within stellar populations (however, see Potter et al. 2012b). The impact of magnetic fields nonetheless may have important consequences on stellar populations and stellar-end products, for example, considering progenitors of magnetars (e.g., Schneider et al. 2019).
We demonstrate that magnetic braking by a fossil field leads to efficient spin-down. For example, an initial equatorial field of 3 kG strength at solar metallicity is sufficient in most models to decrease an initial surface equatorial rotational velocity of 300 km s −1 below 50 km s −1 within the early stages of the main sequence evolution (e.g., Figure 8). For a given magnetic field strength, the spin-down of high-mass stars (> 10 M ) is further aided by mass-loss, whereas the spin-down of lower-mass stars in our grid (< 5 M ) is identifiable due to the long nuclear timescale. The intermediate-mass range (5-10 M ) has the least efficient spin-down over the main sequence evolution.
The "magnetic population" is thus far only identified within the Galaxy by spectropolarimetry and it is unknown what fraction of massive stars possesses strong, surface magnetic fields in extragalactic environments. Generally, the spin-down of the stellar surface for a given magnetic field strength is the most rapid at high metallicity (due to stronger winds), whereas the measurable surface nitrogen abundances are more impacted at lower metallicity (as chemical mixing effects are more pronounced). In the Large Magellanic Cloud, about 20% of stars follow an anomalous pattern on the Hunter diagram, which can be covered with magnetic stellar evolution models. We identify the existence of a range of initial magnetic field strengths (the exact values depending on metallicity, mixing schemes, etc.) that allow for producing slowly-rotating nitrogen enriched Group 2 stars. The lower limit is constrained by a field strength that is needed to brake rotation and produce slow rotators. The upper limit is constrained by a field strength that is needed to allow for rotational mixing and still produce nitrogen enrichment. The range of possible field strengths for the INT models is much narrower than for the SURF models, however, it is compatible with typically measured values. In the LMC and SMC almost all (except some of the lowest mass) models lead to a solution. Contrary, we find that in the Galaxy the formation of Group 2 stars may essentially be prevented for initial masses from 6 to 23 M in the INT/Mix1 scheme ( Figure E14).
Overall, we find significant differences between the braking and mixing schemes. With internal magnetic braking caused by a strong dipolar field, differential rotation cannot develop. With surface magnetic braking caused by a complex magnetic field the physical scenario remains much less clear; the results strongly depend on the chosen assumptions regarding chemical mixing.

ACKNOWLEDGEMENTS
We thank the anonymous referee for providing us with constructive comments. We thank the developers for making their code publicly available. We also acknowledge the open-source code and detailed description of , which helped us with developing our implementation for the SURF models. Furthermore, we greatly appreciate discussions with Henk Spruit, Fabrice Martins, and Philip Dufton. Support for Y.

DATA AVAILABILITY
A full reproduction package is available on Zenodo at https: //doi.org/10.5281/zenodo.7069766, in accordance with the Research Data Management Plan of the Anton Pannekoek Institute for Astronomy at the University of Amsterdam.
The data used in this paper amounts to the order of 1/3 TB. A typical evolutionary model ("history" file in nomenclature) is a few MBs, whereas a typical structure model ("profile" file in nomenclature) is 10 MB. Each evolutionary model has three structure models saved at the ZAMS, mid-MS, and TAMS, respectively. In addition to the evolutionary and structure files, we save and provide the ".mod" files of each run when available. This allows for continuing the computations on the post-main sequence. We also generated isochrones for each sub-grid, these are included in the Zenodo record.
Given the large range of covered parameter space, the output data in this paper is particularly useful for stellar evolution and population synthesis studies, as well as to compare with observational results of even individual stars. However, we emphasise that when interpreting observational results the modelling assumptions and uncertainties should be considered. The MESA EOS is a blend of the OPAL (Rogers & Nayfonov 2002), SCVH (Saumon et al. 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).

APPENDIX B: MASS-LOSS QUENCHING FOR QUADRUPOLE FIELD GEOMETRIES
As a simplifying assumption, we consider that even for more complex magnetic fields (represented by the SURF models), a scaling based on the Alfvén radius of a quadrupole field may be adopted to calculate the mass-loss quenching effect via B (Equations 7-8).
Here, we outline the actual geometrical effects of a quadrupolar magnetic field configuration. Following ud-Doula et al. (2008), the escaping wind fraction is the fraction of the stellar surface covered with open field loops. For a purely dipolar geometry, the region of open field loops corresponds to a polar cap extending from = 1 to = , where is the cosine of the co-latitude , and = cos is the location at which the largest loop connects to the stellar surface (see pink shaded regions in Fig B1).
This assumes that the magnetic field lines retain their pure dipolar geometry, even though in principle the field gets distorted by the outflowing stellar wind. Another approach is to solve for the magnetic field in the magnetosphere using a force-free extrapolation (e.g. Jardine et al. 1999). Assuming that the magnetic field can be expressed as the gradient of a potential Φ such that ∇ 2 Φ( , , ) = 0, the magnetic field can be found by imposing a dipolar (or quadrupolar, see below) boundary condition at the surface and imposing that the magnetic field must become radial at a certain distance from the stellar surface (at the so-called source radius, ). Following the description of Gregory (2011) (their Equation 32 with = 1), in this case the largest field loops (illustrated in the top left panel of B1 with dashed curves) are connected to the surface at: For a quadrupolar configuration, open field loops are located in two polar caps ( cap ) and an equatorial belt ( belt ) (see B1 top panels in blue shading). Thus the total fractional area of open field lines is: For a purely quadrupolar configuration, these footpoints of the largest closed loop can be found by solving for the roots of: and for a quadrupole extrapolated with a source radius: In general, for the same closure radius the escaping wind fraction (surface covered by open field loops) decreases from extrapolated dipole → pure dipole → extrapolated quadrupole → pure quadrupole. We also recall that for a given field strength the closure radius is smaller for a quadrupole field than for a dipole. The lower panels of Fig B1 display the corresponding change in magnetic mass-loss quenching between the dipolar and quadrupolar geometries as a function of closure radius. In a future work we will implement and study magnetic mass-loss quenching for field geometries that deviate from a pure dipole. As shown here, we expect that a quadrupole field will result in a lower escaping wind fraction, i.e., stronger mass-loss quenching for a given magnetic field strength at a given time. However, the time evolution may sensitively impact the results in the evolutionary context.

APPENDIX C: CONVECTIVE EXPULSION AND MAGNETIC SUPPRESSION OF CONVECTIVE LAYERS
"Magneto-convection", the interaction between the magnetic field and turbulent motions in convectively unstable stellar layers, has been extensively studied and recent progress warrants a brief discussion of this topic. Modifying Schwarzschild's criterion for convection, the stability criterion for magneto-convection was first given by Gough & Tayler (1966) and recently revised by MacDonald & Petit (2019), specifically considering radiation pressure to account for sub-surface layers of massive stars, in the form of: where is the thermal expansion coefficient, and are the Alfvén and sound speeds, Γ 1 is the first adiabatic exponent, and is the local pressure. This criterion leads to a critical magnetic field strength, which, if reached, suppresses convection (e.g., Lydon & Sofia 1995;Feiden & Chaboyer 2012), whereas, if it is not reached, the field lines are prevented from reconnecting by strong convective turbulence, which over time leads to diffusing the magnetic flux out of those regions (e.g., Parker 1963). In the one-dimensional stellar evolution models presented in this study, we refrained from using Equation C1. Both MacDonald & Petit (2019) and Jermyn & Cantiello (2020) showed that the key (surface) characteristics during main sequence evolution are not expected to change significantlywhen using Equation C1 -for the stellar mass domain considered in this work, except for the 50-60 M models. We anticipate that, in general, the post-main sequence evolution needs particular attention due to the structural changes which could lead to the development of extended convection zones.
The cores of main sequence massive stars are strongly convective, such that a critical magnetic field strength, which is needed to suppress convection (according to the stability criterion given in Equation C1) may not be reached. In this case, if the magnetic field strength is below its critical strength throughout the core, it is expected that convection will destroy any large scale, organised magnetic flux. This is commonly known as convective expulsion (Parker 1963;Weiss 1966;Spruit 1979;Tao et al. 1998;Schüssler & Knölker 2001). While the convective core should expel any fossil magnetic fields, the development of small-scale magnetic fields via convective dynamos (e.g., Stello et al. 2016;Augustson et al. 2016) may lead to a complex physical picture near the core boundary (Featherstone et al. 2009). It is possible that the fossil magnetic field has a strong impact, for example, by suppressing the radiatively unstable overshooting regions (or even the stellar core, to some extent, see, e.g., Petermann et al. 2015). Thus far, magneto-asteroseismology of two magnetic intermediate-mass stars are consistent with a low value of overshooting (Briquet et al. 2012;Buysschaert et al. 2018). This gives some support of using a modest overshooting, at least, in our INT models, where we assume dipoles embedded deep in the stellar radiative zone. However, this question needs to be studied in more detail in a future work.
A domain of interest is the sub-surface layer of massive stars where inefficient convection (accounting for a few percent of the total luminosity transport) is predicted to arise due to the increase in the opacity of iron, known as iron opacity bumps (Cantiello et al. 2009;Cantiello & Braithwaite 2011;Gräfener et al. 2012). In these radiation-dominated layers, some excess luminosity is transported by convective flux. Using equipartition arguments, Aurière et al. (2007) and Sundqvist et al. (2013) argue that sub-surface convection in radiation-dominated stellar layers can be suppressed by a critical magnetic field when the magnetic pressure reaches the thermal pressure. The exact value of the critical field strength needed to suppress or even inhibit convection remains somewhat uncertaina problem that is mostly explored for the dynamo fields of low-mass stars (Mullan & MacDonald 2001;Feiden & Chaboyer 2012, 2014Feiden 2016;MacDonald & Mullan 2017). Nevertheless, observations clearly show that the Ap phenomenon disappears for stars less massive than about 1.5 M , coinciding with the development of efficient convection (due to the hydrogen opacity bump) in the stellar envelope (Landstreet 1991;Braithwaite & Spruit 2017). We might therefore assume that extended convective zones cannot develop in regions of massive stars where the fossil field (which reaches the critical strength) permeates through. A fossil magnetic field, observable at the stellar surface, should certainly reach sufficient critical strength in these near-surface layers to inhibit convection. All OB stars show evidence for measurable macroturbulent broadening , including magnetic OB stars (Grunhut et al. 2017;Shultz et al. 2018, with the only exception of NGC 1624-2, see Sundqvist et al. 2013), making them a valuable laboratory to further our understanding regarding the origin of macroturbulence. The observed macroturbulence and the putative lack of subsurface convection in magnetic OB stars would support an origin of macroturbulence (at least, in magnetic stars) that is unrelated to sub-surface convection, and could be due to, for example, wave-propagation from the stellar core (Aerts et al. 2009;Bowman et al. 2020).

APPENDIX D: ZAMS AND TAMS STELLAR STRUCTURE MODELS
Here, we describe the stellar structure models of the fiducial model and the corresponding models with initial masses of 20 M in the other braking and mixing schemes at times close to the ZAMS and TAMS. In Figure D1 the INT/Mix1 models are shown. Initially, meridional currents dominate chemical mixing. Given the very efficient magnetic braking in all layers and the solid-body rotation, the angular velocity decreases more than an order of magnitude in all stellar layers from ZAMS to TAMS. Correspondingly, chemical mixing becomes drastically inefficient, leaving a gap between the stellar core and the surface. In Figure D2 the INT/Mix2 models are shown. The results are similar to the INT/Mix1 scheme; however, at the TAMS the effective diffusion coefficient (combining meridional currents and horizontal turbulence) remains to transport chemical elements, although with relatively low values ( ≈ 10 2 cm 2 s −1 ) at the core boundary.
In Figure D3 the SURF/Mix1 models are shown. Similar to the INT/Mix1 scheme, at the ZAMS meridional circulation is the most efficient process driving chemical mixing, while the GSF instability is also present throughout the radiative envelope of the star. At the TAMS, this instability takes over and becomes the dominant process for chemical mixing. Unlike in the INT models, the SURF models show a break in the angular velocity profile, and correspondingly a change in the distribution of specific angular momentum. The stellar core has much more angular momentum (close to its initial value) than in the INT models at the TAMS. In Figure D4 the SURF/Mix2 models are shown. Similarly to the SURF/Mix1 scheme, chemical mixing is initially dominated by meridional circulation before becoming mostly driven by another instability, here by shear mixing. In this scheme, magnetic braking changes the angular momentum profile, and the surface layers have much less angular momentum than the other models. Consequently, the surface rotation of these models at the TAMS is very slow, while again the core maintains a significant fraction of angular momentum. Finally, let us note how the stellar age systematically increases in the TAMS models from the INT/Mix1 (8.56 Myr) to SURF/Mix2 (11.98 Myr) schemes as a consequence of the different braking and mixing assumptions.

E1 Fiducial model on the Hunter diagram
In Figure E1, we follow the evolution of our 20 M fiducial models on the Hunter diagram (Hunter et al. 2008(Hunter et al. , 2009, showing surface nitrogen abundance against the rotational velocity. The ZAMS is at the lower right corner of the diagrams with high surface gravities (colour-coded).
First, we demonstrate how a competition between magnetic braking and rotational mixing alters the evolutionary tracks (see also Paper III). The left panel of Figure E1 shows how the magnetic field strength results in drastic changes of the model predictions. The stronger the magnetic field, the more rapidly rotation brakes, and the less nitrogen can be mixed to the stellar surface. Above a given magnetic field strength, which strongly depends on other stellar and mixing parameters, the braking is efficient enough that nitrogen is not mixed to the surface and the observable nitrogen abundance remains the initial value throughout the main sequence. For the case of the 20 M solar-metallicity model with the INT/Mix1 scheme shown in Figure E1, this threshold occurs roughly at 3 kG.
The right panel of Figure E1 reveals the impact of the initial metallicity on the evolution of the 20 M INT/Mix1 models in the Hunter diagram. The initial equatorial magnetic field strength is chosen to be 1 kG. Consistently with previous findings (see Section 3.2.6), we obtain that chemical mixing is easier to trace in lower metallicity environments where the fractional change between the initial and final abundances are much larger (factor of ∼ 8 at SMC metallicity) than at higher metallicity (factor of ∼ 2 at solar metallicity) for the same fiducial model.

E2 Impact of varying magnetic field strength within the other braking and mixing schemes at solar metallicity
Here, we describe the effect of varying the magnetic field strength for the stellar evolution when accounting for the INT magnetic braking scheme in combination with the Mix2 chemical mixing scheme and the SURF magnetic braking scheme in combination with both the Mix1 and Mix2 chemical mixing schemes at solar metallicity. We display the evolutionary tracks for these models with initial magnetic field strengths of 0.5 and 3 kG in Figures  E2, E3, and E4. The NOMAG/Mix1 models are shown in Figure  9 and the NOMAG/Mix2 models are shown in Figure E2. This is complementary to the description in Section 3.2.5 and Figure 9 for varying the magnetic field strength when using the INT magnetic braking scheme and the Mix1 chemical mixing scheme. Below, we first describe how the stellar evolution is affected by accounting for magnetic fields for each of the SURF/Mix1, INT/Mix2, and SURF/Mix2 cases.

E2.1 SURF/Mix1
As described in Section 3.2 and seen from Figure 9, when the Mix1 rotational mixing scheme is used without including magnetic fields, the stars with initial masses < 30 M remain rapidly rotating (> 200 km s −1 ) until the end of the main sequence. The nitrogen abundance increases smoothly (compared to the initial value, see Table 2) with the initial mass, stretching from log( / ) + 12 = 8.0 up to 8.8. The stars develop mild differential rotation as the main sequence progresses, resulting in cores spinning roughly twice as fast as their envelopes towards the end of the main sequence. When treating magnetic fields using the SURF scheme, angular     momentum is very efficiently transported in the outermost 20% of the mass. In deeper layers of these models efficient shear mixing is invoked by differential rotation (see Section 3.1). As a result, mixing is enhanced in magnetic models treated with the SURF scheme compared to those treated with the INT scheme, which in contrast are rigidly rotating. Nitrogen is mixed out more efficiently, reaching surface abundances up to log( / ) + 12 = 8.7 in the SURF models compared to 8.6 in the INT models of stars with initial equatorial magnetic field strength of 0.5 kG. This can be seen when comparing the middle row of Figure 9 with the top row of E2. In fact, for stars with initial equatorial magnetic field strength of 0.5 kG in the SURF/Mix1 scheme, the surface enrichment of nitrogen is very similar throughout the main-sequence as for the corresponding non-magnetic models. However, the stars spin down somewhat more when magnetic fields are implemented.
In the bottom row of Figure E2, models with initial equatorial magnetic field strength of 3 kG are displayed. From these models it is even more evident that the SURF scheme gives rise to slower spindown than the INT scheme does (cf. the bottom row of Figure 9). While some stars still have a considerable rotational velocity at the end of the main-sequence, magnetic braking is sufficient to produce less enriched stars with nitrogen abundances of log( / ) + 12 = 8.3. When the INT scheme is implemented instead, none of the models with initial magnetic field strength of 3 kG exhibit nitrogenenriched surfaces.

E2.2 INT/Mix2
The Mix2 chemical mixing scheme is an efficient rotational mixing scheme, resulting in quasi-chemically homogeneous main-sequence evolution when magnetic fields are not taken into account (see Section 3.1 and Figure 5). The evolution of the non-magnetic models is displayed in the upper row of panels of Figures E3. We can see that non-magnetic stars are expected to experience modest differential rotation (since a high diffusivity, which we attributed to the magnetic field, is not assumed in these models). All of the non-magnetic models efficiently mix nitrogen out to the surface, reaching values of log( / ) +12 of 8.8-9.5. All but the most massive stellar models remain rapidly rotating throughout most of the main sequence, with surface rotation speeds of ∼ 200 − 350 km s −1 .
In the middle row of Figure E3, we display models initialised with 0.5 kG magnetic field strengths according to the INT scheme. This magnetic field strength is sufficient to brake the rotation of the stars, but only after a significant amount of nitrogen has been mixed out to the surface. As a result, all stars rotate slower than ∼ 50 km s −1 once the TAMS is reached, and stars with initial masses above 7 have surface nitrogen abundances of log( / ) + 12 ∼ 8.9. As visible from the colour-coding in the Kiel diagram, both the most massive and least massive stars spin down early during the main-sequence evolution, while the models in the intermediatemass range of ∼ 5−10 M spin down later during the main-sequence evolution (see Section 3.2.4).
As seen from Figure E3, a modest initial magnetic field strength of 0.5 kG is sufficient to completely alter the evolution in the Hertzsprung-Russell diagram. If magnetic fields are not accounted for, the stars are expected to evolve quasi-chemically homogeneously during the entire main sequence, which is not the case when magnetic fields are taken into account. In the Kiel diagrams, the ratios of the core to envelope spins are shown. When magnetic field effects are accounted for in the INT scheme, the models rotate rigidly throughout the main sequence (cf. Section 3.2.3).
In the case of models with initial magnetic field strength of 3 kG (bottom row of Figure E3), magnetic braking is even more efficient compared to the case of the 0.5 kG models. As a result, rotational mixing can only act briefly before the stars have spun down, leading to slow rotation rates of 50 km s −1 already in the beginning of Figure E2. Same as Figure 9 but for the SURF/Mix1 scheme. Top panels show models with an initial equatorial magnetic field strength of 0.5 kG, whereas the lower panels show models with 3 kG. The NOMAG/Mix1 model is presented in Figure 9. the main sequence (for example when log ∼ 4.0). The maximum surface nitrogen abundances reached are also correspondingly low in comparison to the non-magnetic and weakly magnetic cases, reaching log( / ) + 12 ∼ 8.4.

E2.3 SURF/Mix2
Since stars spin down slower when the SURF scheme is used compared to the INT scheme to implement magnetic fields, the quasichemically homogeneous main sequence evolution that the nonmagnetic models in the Mix2 rotational mixing scheme exhibit is not as easily quenched by the implementation of magnetic fields in the SURF/Mix2 combination. The effect of the weaker magnetic braking is evident when comparing Figures E3 and E4.
As a result of the slower spin-down, even stars with initial magnetic field strength of 3 kG still substantially enrich their surfaces with nitrogen (log( / ) + 12 = 8.7 − 9.2). They also reach rotation rates close to zero at the end of the main sequence. As seen in Figure 5, the combination of SURF and Mix2 results in such efficient mixing that even the surface helium abundance significantly increases during the main sequence.
E3 Impact of varying magnetic field strength at = 0.0064

E3.1 INT/Mix1
At LMC metallicity, the non-magnetic models have several main differences compared to solar metallicity models. First, at lower metallicity the position of the ZAMS is shifted to higher effective temperatures. Second, the relative amount of surface nitrogen enrichment is much higher. For comparison, the 60 M model in the Mix1 chemical mixing scheme at solar metallicity produces an order of magnitude enrichment in its surface nitrogen to hydrogen ratio over its main sequence evolution, whereas the LMC model produces almost a factor of 30 (see Figure E5).
With a 0.5 kG initial equatorial magnetic field strength within the INT/Mix1 scheme, the surface equatorial rotational velocity is reduced in the models due to magnetic braking. This produces slight differences in the evolutionary tracks on the HR and Kiel diagrams. As a consequence of the reduced rotation in the star, chemical mixing is also somewhat less efficient compared to the non-magnetic case.

E3.2 SURF/Mix1
With a 0.5 kG initial equatorial magnetic field strength within the SURF/Mix1 scheme ( Figure E6), the overall reduction in surface equatorial rotational velocities is somewhat less than in the INT/Mix1 case. On the other hand, the most massive models in this configuration tend to achieve lower rotational velocities at the TAMS.
With a stronger, 3 kG initial equatorial magnetic field strength, the numerical simulations show that the position of the models on HR and Kiel diagram are similar to the non-magnetic and 0.5 kG cases, except for the most massive stars, which can reach much lower surface gravities towards the end of their main sequence evolution. As expected, the evolution of surface equatorial rotational velocities is also impacted, showing that a spin-down is achieved in all models.
(However, see Section 3.2.4 for the discussion on models in the mass range of 8-15 M .) Consequently, the surface nitrogen enrichment is more modest compared to the non-magnetic case but still notable in these models.

E3.3 INT/Mix2
When chemical mixing is assumed in the Mix2 scheme, the nonmagnetic models at LMC metallicity evolve bluewards, quasichemically homogeneously for their entire main sequence evolu- Figure E4. Same as Figure 9 but for the SURF/Mix2 scheme. See top panels of Figure E3 for the NOMAG/Mix2 models. tion. Within this configuration the [ / ] ratio reaches extremely high values ( Figure E7). This is in part due to the efficient mixing of nitrogen, and in part due to losing a significant amount of hydrogen from the stellar surface via winds. In this case, the low surface rotational velocities at the TAMS are achieved due to Wolf-Rayet type mass loss taking away angular momentum from the star.
When a 0.5 kG initial equatorial magnetic field strength is considered within the INT/Mix2 scheme the models return to a redward evolution on the main sequence following a rapid blueward evolution. This is because the spin-down due to magnetic braking reduces the overall efficiency of chemical mixing and prevents a quasi-chemically homogeneous evolution for the entire main sequence. Note that in this case the TAMS surface rotational velocities are actually higher than in the non-magnetic case as magnetic braking weakens over time and stellar winds are less efficient compared to the non-magnetic case.
With a 3 kG initial equatorial magnetic field strength, the main differences result from a stronger magnetic braking. This almost immediately prevents a blueward evolution, producing slowly-rotating tracks in which the surface nitrogen enrichment is significantly reduced compared to the non-magnetic case.

E3.4 SURF/Mix2
In the SURF/Mix2 scheme both set of models with 0.5 and 3 kG initial equatorial magnetic field strength undergo a rapid initial blue-ward evolution as in the INT/Mix2 case ( Figure E8). With weaker initial magnetic fields, the models above10 M turn to a short redward evolution before spinning down and reaching their TAMS again at hotter effective temperatures. With stronger initial magnetic fields, the spin-down leads to preventing the models from turning back to a blueward evolution by the end of their main sequence. Similarly to the SURF/Mix1 case, the most massive models with 3 kG initial equatorial magnetic field strength may reach low surface gravities. In this configuration, it is also possible that the most massive models produce a final [ / ] ratio that is lower than that of lower mass models (see also Figure 11). However, given the very efficient Mix2 scheme, the predicted surface nitrogen enrichment is still expected to be observable as it is over an order of magnitude compared to the baseline. E4 Impact of varying magnetic field strength at = 0.0026

E4.1 INT/Mix1
At SMC metallicity, the ZAMS position of the non-magnetic models is further shifted to higher effective temperatures compared to solar and LMC metallicity models ( Figure E9). For the 50 and 60 M models, we note major differences in comparison to models at higher metallicity. The highest mass non-magnetic models at SMC metallicity are able to undergo a blueward evolutionary phase for some fraction of their main sequence. This results from the relatively Figure E5. Same as Figure 9 but for the NOMAG/Mix1 scheme (top) and INT/Mix1 scheme with an initial equatorial magnetic field strength of 0.5 kG (lower panel) at = 0.0064. The INT/Mix1 model with initial equatorial magnetic field strength of 3 kG is presented in Figure 10. rapid rotation of the models, the efficient core-envelope coupling (assumed via hydrodynamic processes in a diffusive scheme), and the effects of stellar winds. These models can enhance their surface [ / ] ratio by about a factor of 50 compared to the initial baseline value.
When a 0.5 kG initial equatorial magnetic field strength is considered, the effects of the spin-down are more pronounced than in higher metallicity models. We see that the blueward evolution of the most massive models is prevented. The impact on the nitrogen abundances is also notable as less enrichment is predicted compared to non-magnetic models. However, we should also note that in contrast to the most massive models in the INT/Mix1 scheme at higher metallicities, the rotational velocities of these models at SMC metallicity is actually less impacted. This is because stellar winds are less efficient at lower metallicity. Thus the combination of magnetic braking by a "weak" field (initially 0.5 kG) and weaker stellar winds allows for maintaining higher surface rotational velocities.

E4.2 SURF/Mix1
In the SURF magnetic braking scheme we notice small differences compared to the INT scheme when keeping chemical mixing the same in the Mix1 scheme ( Figure E10). For the same initial mass the SURF models tend to maintain a somewhat higher rotational velocity than the INT models. As expected, when increasing the value of the initial equatorial magnetic field strength from 0.5 to 3 kG, the spin-down is stronger, and consequently lower surface rotational velocities can be reached in conjunction with a reduction in the surface nitrogen enrichment. Nonetheless, the most massive models in this configuration with a 3 kG initial equatorial magnetic field strength are able to produce slowly-spinning, nitrogen-enriched Group 2 stars.

E4.3 INT/Mix2
In the NOMAG/Mix2 at SMC metallicity, the most extended blueward evolution is produced, reaching 80 kK for the most massive stellar models ( Figure E11). The behaviour of these models is qualitatively similar to those of solar and LMC metallicity. The predicted increase in [ / ] is first a factor of about 50 related to mixing nitrogen, and then further almost a factor of 10 related to losing hydrogen from the stellar surface. Due to the blueward evolution of the models the surface gravities are expected to remain high.
In the INT/Mix2 case, a 0.5 and 3 kG initial equatorial magnetic field strength can prevent the extensive blueward evolution of the models. However, in contrast to some previous cases, the most massive models do not reach low surface gravities. Although for the non-magnetic case, all models are predicted to show a very large increase in their surface nitrogen abundance, the magnetic models show a more moderate enrichment. In this particular setup Figure E6. Same as Figure E5 but for the SURF/Mix1 scheme. Top panels show models with an initial equatorial magnetic field strength of 0.5 kG, whereas the lower panels show models with 3 kG. The NOMAG/Mix1 model is presented in Figure E5.
we also see that models less massive than initially 6 M tend to produce much lower [ / ] ratios, whereas models with higher initial masses are closely grouped together.

E4.4 SURF/Mix2
In the SURF/Mix2 scheme at SMC metallicity the 0.5 kG initial magnetic field strength is not quite strong enough to prevent the quasi-chemically homogeneous evolution of the models ( Figure  E12). On the other hand, a 3 kG initial equatorial magnetic field strength leads to a more pronounced redward evolution of the most massive models after their spin-down. In this case, we again observe that the most massive models can reach low surface gravities along with low surface rotational velocities. In all cases, the predicted enrichment in [ / ] remains over an order of magnitude compared to the baseline.

APPENDIX F: CUTOFF MAGNETIC FIELD STRENGTH IN THE SMC AND SOLAR NEIGHBOURHOOD
As in Section 4.2, we can quantify the cutoff magnetic field strengths in other metallicity environments needed to achieve slow-rotating, nitrogen-enriched stars. First, let us turn our attention to models at the SMC, shown in Figure E13. The four schemes predict essentially constant values for max,N for masses higher than about 10 M . The INT/Mix1 scheme is the most limiting scenario from the four schemes. The other schemes predict a wider range of magnetic field strengths that allow for producing Group 2 stars. Since nitrogen abundances and rotation rates are measurable in the SMC, the identified Group 2 stars (Dufton et al. 2020) could be explained with a limited range of magnetic field strengths. Typically, the magnetic field strength decreases roughly by an order of magnitude on the main sequence (cf. Equation 2). Thus depending on the age and evolutionary state of observed Group 2 stars, their observable fields will be weaker than the ZAMS field strengths. A sort of minimum range of their putative magnetic field strength could evolve from the range of initially 250 G to 7.5 kG (equatorial strength) at the ZAMS to 25 G to 750 G at the TAMS if the INT/Mix1 scheme is considered. For massive star models, the SURF/Mix1 and SURF/Mix2 schemes allow for an upper limit of initially 50 kG, which we expect to weaken to a 5 kG equatorial field strength by the TAMS. Figure E14 shows models with an initial metallicity of = 0.014 (representative of the Solar neighbourhood). Consistently with previous findings, the INT/Mix1 is the most restrictive to produce Group 2 stars. In fact, in this case there are more forbidden regions than in models with lower metallicities. For example, the INT/Mix1 scheme essentially does not allow for producing Group 2 stars in the 6 -23 M range. In contrast, the SURF/Mix2 scheme allows for strong magnetic fields to brake the rotation and still produce nitrogen excess at the surface. Since high-resolution spectropolarimetric observations in the Galaxy have allowed for measuring the magnetic field strengths of massive stars, these model predictions allow for more direct comparison with observations than in the Magellanic Clouds. While the nitrogen abundance and (projected) rotational velocity are also measurable, a comprehensive study still needs to assess these three quantities jointly in comparison to magnetic stellar evolution models (however, see Morel et al. 2008;Martins et al. 2012Martins et al. , 2015Aerts et al. 2014). Our single-star models provide testable predictions for these parameters, and such an in-vestigation could help disentangle between the uncertain braking and mixing schemes. In particular, nitrogen-enriched slow rotators with strong magnetic fields favour efficient mixing. This paper has been typeset from a T E X/L A T E X file prepared by the author. Figure E14. Same as Figure E13 but for Solar metallicity models.