© The Authors 2025

1. Introduction

The discovery of luminous quasars hosting billion-solar-mass supermassive black holes (SMBHs) at redshifts z ≳ 7 poses significant theoretical challenges to our understanding of black hole formation and growth in the early universe (Mortlock et al. 2011; Wu et al. 2015; Bañados et al. 2018). The existence of these massive objects within the first billion years necessitates rapid formation pathways, motivating models invoking heavy seed black holes originating from the direct collapse of primordial gas into supermassive stars (SMSs) in the mass range of 10³–10⁵ M_⊙ (Latif et al. 2016; Smidt et al. 2018; Inayoshi et al. 2020; Latif et al. 2022; Patrick et al. 2023). Complementary evidence from recent JWST observations further complicates this picture, revealing galaxies at z > 8 with anomalously high nitrogen-to-oxygen ratios (N/O), significantly exceeding standard enrichment scenarios (Cameron et al. 2023; Nakajima et al. 2023; Isobe et al. 2023; Larson et al. 2023). These observations suggest the early universe hosted exotic stellar populations capable of rapid chemical enrichment (Bunker et al. 2023; Marques-Chaves et al. 2024; Liu et al. 2025).

Several mechanisms have been proposed to explain these chemical anomalies, including enrichment by Wolf–Rayet starbursts, asymptotic giant branch (AGB) stars, and extremely massive stars (Vink 2023; Nandal et al. 2024c). Wolf–Rayet stars and AGB stars, while prolific nitrogen producers, typically require metallicities or evolutionary timescales that are incompatible with observed high N/O values at such early epochs (Cameron et al. 2023). Alternatively, SMSs are compelling candidates due to their potential to synthesize and rapidly disperse large amounts of primary nitrogen, produced through efficient CNO cycling in their cores, into their surrounding environments via strong mass-loss events or explosive endpoints (Nagele & Umeda 2023; Nandal et al. 2025a).

Hydrodynamical simulations have demonstrated viable conditions for SMS formation possibly via strong external Lyman–Werner radiation (Latif et al. 2014; Schauer et al. 2017; Regan et al. 2017) or dynamically heated environments (Regan et al. 2017; Wise et al. 2019; Regan et al. 2020b), or through the suppression of fragmentation by baryon-dark matter streaming motions (Schauer et al. 2017; Liu et al. 2024; Ishiyama & Hirano 2025). More recent simulations also show that rapid mass growth via cold accretion flows can lead to SMS formation (Latif et al. 2022; Kiyuna et al. 2023, 2024). Such conditions permit the formation of very massive objects and subsequently massive seed black holes, directly addressing the SMBH timing problem (Woods et al. 2017; Whalen et al. 2020). In instances where the host dark matter halo can grow to sufficient masses, the gas inflow rate into the central regions of the halo can exceed 10⁻²–10⁻¹ M_⊙ yr⁻¹, allowing SMSs to form in such halos. An alternative pathway that facilitates such high accretion rates is the Population III.1 (Pop III.1) model, which invokes heating from weakly interacting massive particle (WIMP) annihilation (Spolyar et al. 2008; Tan et al. 2024). In this scenario, the slow collapse of gas within an isolated pristine mini-halo can adiabatically contract the co-located dark matter, thus boosting its density to a level where WIMP annihilation becomes a dominant energy source for the forming protostar (Banik et al. 2019). This energy injection supports the protostar in an inflated ‘swollen’ state, maintaining a cool photosphere (∼10⁴ K) that emits few ionizing photons (Nandal et al. 2025a). The resulting suppression of radiative feedback allows the star to bypass the mass limit typical of standard Pop III stars and grow efficiently, potentially reaching supermassive scales of ∼10⁵ M_⊙ by accreting the bulk of its host mini-halo’s baryonic content (Singh et al. 2023).

In some cases (where high accretion rates can be maintained; Latif et al. 2021), protostellar growth to ∼10⁵ M_⊙ can be achieved before the general-relativistic instability triggers collapse to a massive black hole seed (Banik et al. 2019; Wise et al. 2019; Regan et al. 2020a; Regan 2023; Singh et al. 2023; Herrington et al. 2023). However, current stellar evolution models for SMSs assume no rotation, and thus neglect centrifugal support and rotationally induced chemical mixing. State-of-the-art stellar evolution models also omit rotationally enhanced mass loss and angular momentum transport that determines black hole spin, precluding theoretical predictions of remnant spin distributions (Woods et al. 2020; Nandal et al. 2023).

Rotation significantly influences massive star evolution by enhancing internal mixing through meridional circulation and shear instabilities, enlarging convective cores, and extending stellar lifetimes (Maeder 1997; Heger et al. 2000; Ekström et al. 2012; Nandal et al. 2024b). In standard massive stars (∼10–120 M_⊙), rotation can induce chemically homogeneous evolution, drastically altering their evolutionary pathways and final fates (Yoon & Langer 2005; Yoon et al. 2012). Additionally, rotational mixing affects stellar surface abundances, enriching atmospheres with processed elements such as helium and nitrogen much earlier than in non-rotating stars (Brott et al. 2011; Limongi & Chieffi 2018; Murphy et al. 2021; Tsiatsiou et al. 2024; Nandal et al. 2024d). Rotation also enhances mass-loss rates through mechanical winds, particularly near critical rotation thresholds, thereby significantly influencing stellar mass evolution and the resulting remnant masses (Langer 1998; Vink 2018).

Near the Eddington limit, the interplay between radiation pressure and rotation results in the ΩΓ limit, a threshold where effective gravitational acceleration at the stellar surface nearly vanishes (Langer 1998; Maeder & Meynet 2000). Stars approaching this limit experience drastically enhanced mass loss, even at rotation rates below classical break-up velocities (Glatzel 1998; Owocki et al. 2004). This rotational-radiative instability is particularly relevant for SMSs, given their inherently high luminosities near the Eddington limit, potentially imposing strict caps on their rotation and significantly influencing their evolution and ultimate fates (Owocki et al. 2004; Sanyal et al. 2015).

Rotation critically impacts the angular momentum content and resulting spins of black hole remnants from massive stars (Heger et al. 2000). Intermediate-mass black holes (IMBHs), including those from SMS collapse, inherit angular momentum profiles established during stellar evolution (Woosley & Heger 2006; Fryer & Heger 2011). Black hole spins strongly affect accretion disk dynamics, jet production, and gravitational-wave emission during mergers, thereby influencing their observational signatures and theoretical interpretations (Berti & Volonteri 2008; Barausse 2012). Understanding the rotational evolution of progenitor stars thus provides crucial constraints on IMBH spins, with direct implications for SMBH formation scenarios through hierarchical merging (Sesana et al. 2009).

To date, no models of SMSs above 1000 M_⊙ with rotation have been created, leaving substantial uncertainties in the predictions of the evolution of these SMSs and their ultimate fates. Motivated by observational puzzles and theoretical imperatives, we present in this paper the first comprehensive set of rotating evolutionary models for Population III SMSs with masses of 1000–10 000 M_⊙. Our models incorporate detailed treatments of rotation-induced mixing, angular momentum transport, and the effects of the ΩΓ limit.

In Section 2 we describe our numerical methods and input physics, including the implementation of rotation in the stellar evolution code GENEC. Section 3 presents detailed evolutionary results, highlighting the influence of rotation on stellar lifetimes, internal structures, and mass loss. In Section 4 we compare our rotating SMS models to previous non-rotating calculations and discuss their implications for early chemical enrichment and SMBH seed formation scenarios. Finally, in Section 5 we summarize our key findings and outline promising directions for future research.

2. Methods

We present the first rotating models of supermassive (1000-10 000 M_⊙) Pop III stars and investigate the effects of the ΩΓ limit on their evolution. These models use the same numerical set-up and physics prescriptions as in prior non-rotating studies (Nandal et al. 2024a), differing only by the inclusion of rotation and thereby isolating its effects. All the models were computed using the Geneva Stellar Evolution Code (GENEC; Eggenberger et al. 2008; Ekström et al. 2012; Nandal et al. 2024b).

2.1. Initial ZAMS structures

For stellar masses exceeding about 1000 M_⊙, the notion of a classical ZAMS is not well defined. We therefore began by constructing non-rotating, accreting pre-MS Pop III models using the tracks of Nandal et al. (2025a). These tracks begin at 2 M_⊙ and, under constant accretion rates, grow until the target mass (between 1000 M_⊙ and 10 000 M_⊙) is attained. While accretion in nature is highly variable (Regan et al. 2020a), the final ZAMS structure is primarily determined by the star’s total mass and composition, as it reaches thermal equilibrium during the final contraction phase, largely erasing the memory of its prior accretion history (Nandal et al. 2023). Despite remarkable successes in the 3D hydrodynamic simulations that account for angular momentum transport in the pre-MS (Hirano & Bromm 2018; Kimura et al. 2023), a self-consistent 1D prescription for stellar evolution codes does not yet exist in the literature. Therefore, our pre-MS models do not include rotation. In addition, this allowed us to generate a consistent set of non-rotating ZAMS models to serve as a baseline for our study.

Rotation was then introduced with surface velocities of 0.01, 0.05, and 0.10 of the critical Keplerian velocity (v/v_crit). These relatively low rates were chosen specifically to test the regime where the ΩΓ limit is predicted to become important for these massive stars (Maeder & Meynet 2000; Haemmerlé 2021) and to quantify its effect from the onset of the main sequence. The construction of these initial rotating models is detailed below.

The initial structure for each rotating model was generated using a dedicated numerical procedure with three key steps:

1. Contraction to the ZAMS: Once a non-rotating pre-MS model reaches its final target mass, accretion is terminated. The model contracts on a Kelvin–Helmholtz timescale until it reaches its highest effective temperature, at which point its structure is in near-thermal equilibrium. The resulting composition was X = 0.7516; Y = 0.2484; Z = 0. The physical parameters of this non-rotating ZAMS model, including the profiles of pressure, temperature, density, and chemical composition are then extracted.

2. Incorporation of Rotation: The extracted non-rotating structure served as the input for a new GENEC simulation via a dedicated FORTRAN routine. At the start of this run, angular momentum was injected to produce the desired surface rotation rate (e.g. 0.1 v/v_crit). The mapping from the non-rotating structure to the rotating one is handled by a dedicated FORTRAN wrapper that passes the full radial grids of P; T; ρ; and X_i to GENEC. The code then self-consistently evolves the transport of angular momentum and chemical species.

3. Relaxation to a Rotating Equilibrium: After rotation is introduced, the model evolves for several thousand years. This timescale is negligible compared to the main-sequence lifetime but is sufficient for the star to dynamically and thermally relax to a stable, rotating equilibrium. We define this converged state as the rotating ZAMS. The model settles on the blue side of the HR diagram at its maximum T_eff, which we adopt as the ZAMS reference point. This procedure ensures our investigation begins from a physically consistent starting point.

The main sequence is defined to begin when the central hydrogen mass fraction, X_c, has decreased by 10⁻³ from its initial value, and to end when X_c falls below 10⁻³. These self-consistent initial models set the stage for our rotational prescription, described next.

2.2. Treatment of rotation

We adopted the shellular rotation framework of Zahn (1992), in which strong horizontal turbulence enforces nearly uniform angular velocity on isobaric surfaces. The transport of chemical species is described by a purely diffusive equation:

$$\begin{array}{c}\hfill \rho \frac{\partial X_i}{\partial t}{|}_{M_r}=\frac{1}{r^2}\frac{\partial }{\partial r}\left(\rho r^2{D}_{\mathrm{chem}}\frac{\partial X_i}{\partial r}\right).\end{array}$$(1)

Here X_i is the mass fraction of species i,

$$\begin{array}{c}\hfill D_{\mathrm{chem}}=D_{\mathrm{shear}}+D_{\mathrm{eff}},\end{array}$$(2)

and

$$\begin{array}{c}\hfill D_{\mathrm{eff}}=\frac{1}{30}\frac{[r{U}_2(r){]}^2}{D_{\mathrm{h}}}.\end{array}$$(3)

We took D_shear from Maeder (1997) and the horizontal turbulence D_h from Zahn (1992). Both prescriptions are standard in Geneva-type models of massive stars up to 120 M_⊙.

The angular momentum distribution is evolved through the Lagrangian form

$$\begin{array}{c}\hfill \rho \frac{\partial }{\partial t}{\left(r^2\mathrm{\Omega }\right)}_{M_r}=\frac{1}{5r^2}\frac{\partial }{\partial r}[\rho r^4\mathrm{\Omega }{U}_2(r)]+\frac{1}{r^2}\frac{\partial }{\partial r}[\rho D_{\mathrm{ang}}{r}^4\frac{\partial \mathrm{\Omega }}{\partial r}],\end{array}$$(4)

where Ω is the (shellular) angular velocity, U₂(r) denotes the radial component of the meridional circulation, and D_ang = D_shear for these non-magnetic models.

Additional physics is captured by the effective diffusion coefficient D_eff, which accounts for the mixing induced by the interplay of meridional circulation and horizontal turbulence. Our prescription for D_shear models the destabilizing effect of differential rotation countered by the stabilizing influence of thermal and compositional gradients, leading to efficient mixing in the radiative zones. Furthermore, the assumption of strong horizontal turbulence smooths out latitudinal variations, allowing a one-dimensional (radial) treatment of the mixing. These rotational transport mechanisms critically affect both the chemical stratification and angular momentum distribution, thereby influencing the evolutionary paths of our supermassive star models.

Line-driven winds at Z = 0 were included using the H and He resonance line prescription of Kudritzki (2002) and Ekström et al. (2008) for log(L/L_⊙) > 6. The formalism retains the canonical Ṁ ∝ (Z/Z_⊙) 0.5 scaling and channels photon momentum through the strongest hydrogen and helium transitions. This yields rates far below those of metal-rich counterparts yet provides a physically motivated sink of surface angular momentum. This results in extremely low mass-loss rates, where even the most massive 9000 M_⊙ model only loses ∼2 M_⊙ by the end of the main sequence. Pulsation-driven mass loss is poorly constrained at zero metallicity and is therefore excluded. The models thus evolve at almost constant mass until the surface reaches the ΩΓ limit and rotational–radiative coupling drives outflows.

2.3. The ΩΓ limit

Massive stars approaching the Eddington limit experience a reduction of their effective gravity due to the combined effects of centrifugal and radiative forces (Langer 1998; Glatzel 1998; Maeder & Meynet 2000). In our shellular rotation framework the total effective gravity is written as

$$\begin{array}{c}\hfill g_{\mathrm{tot}}=g_{\mathrm{eff}}(1-\mathrm{\Gamma }),\end{array}$$(5)

where g_eff = g_grav + g_rot and g_grav = GM/r² is the gravitational acceleration (with gravitational constant G, stellar mass M, and local radius r), and g_rot is the centrifugal acceleration. Here Γ represents the local Eddington factor, which compares the actual radiative flux to its maximum allowed value and may vary with latitude due to gravity darkening.

Setting g_tot = 0 leads to two distinct critical velocities. The first critical velocity is obtained by neglecting radiative forces, so that g_eff = 0, which yields the classical break-up condition,

$$\begin{array}{c}\hfill v_{\mathrm{crit},1}=\sqrt{\frac{2GM}{3R_{\mathrm{pb}}}},\end{array}$$(6)

where R_pb is the polar radius at break-up. In the presence of significant radiative acceleration, a second solution arises by imposing g_eff (1 − Γ) = 0, i.e. when Γ = 1. This yields a lower critical velocity that accounts for both rotational and radiative effects. Following Maeder & Meynet (2000), the second critical velocity is

$$\begin{array}{cc}\hfill v_{\mathrm{crit},2}^2& ={\mathrm{\Omega }}^2{R}_{\mathrm{e}}^2(\omega )\hfill \\ \hfill & =\frac{81}{16}\frac{[1-{\mathrm{\Gamma }}_{max}]GM}{V\prime (\omega )R_{\mathrm{eb}}^3(\omega )}\hfill \\ \hfill & =\frac{9}{4}{v}_{\mathrm{crit},1}^2[1-{\mathrm{\Gamma }}_{max}]\frac{R_{\mathrm{eb}}^2(\omega )}{R_{\mathrm{e}}^2(\omega )},\hfill \end{array}$$(7)

where Ω is the angular velocity; R_e(ω) is the equatorial radius for a dimensionless rotation parameter ω; R_eb(ω) is the equatorial radius at break-up; Γ_max is the maximum Eddington factor over the stellar surface; and V′(ω) is a dimensionless function of ω (see Maeder & Meynet 2000 for further details). In the ΩΓ limit, defined by the condition Γ = 1, even moderate rotation drives the star to break up, leading to dramatically enhanced mass loss. This second solution, v_crit, 2, becomes physically relevant for stars with Γ_max ≳ 0.639, a regime typical of our supermassive star models.

In our models, when the surface equatorial velocity, v_eq, approaches this second critical velocity, v_crit, 2, the star enters a phase of intense mechanical mass loss. This process is not treated with a continuous wind prescription but is implemented as a numerical enforcement of the physical limit. The code continuously monitors the ratio v_eq/v_crit, 2. If this ratio reaches unity, a mass-loss routine is activated that removes the necessary amount of mass (ΔM) and angular momentum from the surface layer in that timestep to ensure the star remains just below the critical limit in the subsequent step. The mass-loss rate is thus dynamically calculated by the code to enforce the condition

$$\begin{array}{c}\hfill \dot{M}=\{\begin{array}{cc}{\dot{M}}_{\mathrm{needed}}\hfill & \text{if}{v}_{\mathrm{eq}}\ge v_{\mathrm{crit},2}\hfill \\ 0\hfill & \text{if}{v}_{\mathrm{eq}}<v_{\mathrm{crit},2}\hfill \end{array},\end{array}$$(8)

where Ṁ_needed is the rate required to maintain the star at the stability boundary. This self-regulating mechanism effectively caps the surface rotation of the star.

3. Results

We show the evolution of the three SMSs at v/v_crit = 0.01, 0.05 and 0.1 in the HR diagrams in panels a, b, and c of Figure 1 respectively. Total lifetimes (column 3), times spent as blue supergiants (column 4) and red supergiants (column 5) are shown in Table 1 along with convective core mass fractions (columns 6 and 10), equatorial velocities (columns 7 and 11), total angular momenta (columns 8 and 12) and Eddington factors (columns 9 and 13).

Fig. 1. SMS evolution on the HR diagram. The colourbar indicates central H mass fractions, Xc, and the shaded red regions show the isoradii lines: the lightest shade corresponds to log (R/R⊙) = 1 and the darkest shade to log (R/R⊙) = 4. The white outline around the stellar tracks correspond to the evolutionary stage where the ΩΓ-limit is encountered. Top row: 1000–10 000 M⊙ stars at initial rotation rates v/vcrit = 0.01 (a), 0.05 (b), and 0.10 (c). Bottom row: comparison of rotation rates at a given mass, 1000 M⊙ (d), 5000 M⊙ (e), and 9000 M⊙ (f).

3.1. Evolution in the HR diagram

All three sets of models begin their lifetimes as hot and compact blue supergiants with stellar stellar radii of 16–55 R_⊙ and an effective temperature around 88 000 K. Following a short contraction phase that lasts a few thousand years, all models reach a central temperature of log (T_eff) = 8.15, marking the start of core hydrogen burning. All models possess large convective cores that make up about 95% of the total star in mass. As hydrogen continues to fuse in the convective core, the total energy generated is initially dominated by the pp chain. However, the central temperature is high enough for the 3-α reaction to produce ¹²C and once its central abundance exceeds a mass fraction of 10⁻¹⁰, the CNO cycle is initiated. The energy generation from CNO burning contributes significantly to the total energy production of the star, and this energy is transported from the centre to the surface. A fraction of this energy goes into in expanding the envelope, thereby increasing the stellar radius. As the models continue to deplete their central hydrogen reservoir, their stellar radius increases as they steadily begin their transition into red supergiant protostars with stellar radii varying from 120 to 10 000 R_⊙.

The expansion of the stars naturally affects their equatorial velocities throughout their evolution. For instance, the equatorial velocities range between 26 and 44 km s⁻¹ at the start of core hydrogen burning, making these particular models slow rotators in spite of their compact size. This is simply due to the choice of initial rotation, which is based on the ΩΓ limit, and is explored in a later section. These equatorial velocities later fall to nearly zero, especially in case of more massive models (5000–10 000 M_⊙) as the central hydrogen is depleted. An interesting feature to note here is the position of the tracks in the HR diagram at the end of core hydrogen burning. We find that as SMS mass increases, core hydrogen burning ends in a more red-ward position in the HR diagram. This result is in agreement with standard rotating 9–120 M_⊙ star models. The main-sequence lifetime of these stars ranges from 1.0 Myr (10 000 M_⊙) – 1.6 Myr (1000 M_⊙), as shown in column 2 of Table 1. Although the mass difference between the two models is one order of magnitude (1000%), the difference in age is only 60%. This is a due to a shift in the luminosity-mass relationship, where the scaling shifts to L ∝ M. In addition, all models spend the majority of their lifetimes (90–100%) on the blue side of the HR diagram (log (T_eff) > 4), as seen in columns 4 and 5 of Table 1.

As shown in panel b of Figure 1, the 0.05 v/v_crit SMS exhibits evolutionary trends in nuclear reaction rates, effective temperatures, and convective core properties that are quite similar to those of the 0.01 v/v_crit model. The slight but important differences arise in the total main-sequence lifetimes, as they range from 1 Myr (10 000 M_⊙) – 1.74 Myr (1000 M_⊙), as shown in column 2 of Table 1. This increase in lifetime on the lower-mass end, for instance at 1000 M_⊙, is around 7% larger than the same star at the lowest rotation rate. This effect is primarily due to stronger elemental mixing (in this case hydrogen) during the main sequence made possible by the increased rotation rate. Hydrogen from the radiative envelope is dredged down to the convective core by rotational mixing facilitated by meridional currents, which partially resupplies the core and extends its size. This can be seen in columns 6 and 10 of Table 1, where the convective core mass fraction, M_cc, is larger for the 0.05 v/v_crit models than for the corresponding 0.01 v/v_crit models, especially by the end of core hydrogen burning. Because of the higher rotation rate, the equatorial velocities (columns 7 and 11 of Table 1) are naturally higher for 0.05 v/v_crit models at the onset of core hydrogen burning, but this trend is not seen at the end of core hydrogen burning in the 5000–10 000 M_⊙ models. These stars also approach the Hayashi limit and their radii, and therefore their equatorial velocities, are identical to those of the 0.01 v/v_crit models. The total angular momentum L_tot reservoirs (columns 8 and 12 of Table 1) of these models are also higher throughout their evolution. A point of similarity between the two sets of models here is the Eddington factor (Γ_edd), shown in columns 9 and 13 of Table 1. All models, independent of rotation rates, exhibit very high Γ_edd, which makes them susceptible to the ΩΓ limit (see Section 3.4).

The general evolutionary trends of the 0.1 v/v_crit model shown in panel c of Figure 1 are again similar to those at 0.01 and 0.05 v/v_crit, with a few additional insights. The extent of chemical mixing in these models is highest among the three rotation rates and produces the largest M_cc. This in turn extends their lifetimes by another few percent when compared to the 0.05 v/v_crit models. Stronger mixing also boosts energy production in the stars. As this energy propagates through the stellar interior it expands the envelope to even greater radii than those of the other two sets of models, ensuring that stars at these highest rotation rates spend the longest times in the red (column 5 of Table 1) and subsequently end the main sequence on the red side of HR diagram. Although these models have the highest initial equatorial velocities and total angular momentum reservoirs, these velocities drop to zero once they approach the Hayashi line and become red supergiants.

We compare the evolution of 1000 M_⊙, 5000 M_⊙, and 9000 M_⊙ stars at all three rotation rates in panels d, e, and f of Figure 1. The effect of rotational mixing is evident in the 1000 M_⊙ HR tracks when they begin to diverge after beginning the main sequence at the same place in the HR diagram. The 0.01 v/v_crit star exits the main sequence with log (L/L_⊙) = 7.51 and log (T_eff = 4.6. As chemical mixing is enhanced by higher rotation rates, the 0.05 v/v_crit star becomes slightly more luminous, log (L/L_⊙) = 7.53 and cooler, log (T_eff) = 4.5. This effect is enhanced in the 0.1 v/v_crit model, where log (L/L_⊙) = 7.56 and log (T_eff) = 3.8. At log (L/L_⊙) = 7.525, the 0.1 v/v_crit model fluctuates between the red and blue sides of the HR diagram because Γ_edd exceeds 0.9, causing the star to encounter the ΩΓ limit. The consequences of this limit will be discussed later.

The effects of different rotation rate are also evident in the 5000 M_⊙ models in panel e of Figure 1, as the luminosities vary by ∼ log (L/L_⊙) = 0.075. The effective temperature of the 0.05 v/v_crit SMS is much higher than for the other two models, as it ends its core hydrogen burning stage at log (T_eff) = 4.4. Generally, this effect can be attributed to stronger mixing in the convective core induced by rotation and this fact is true when comparing the 0.01 v/v_crit and 0.1 v/v_crit models. The 0.1 v/v_crit star is indeed redder because of stronger rotational mixing, which allows more efficient transport of hydrogen into the core and builds up the helium core at a higher rate. This raises the luminosity even higher and the model approaches the Eddington luminosity sooner. But 0.05 v/v_crit being less red than 0.01 v/v_crit could be attributed to the core properties and temperature sensitivity of the CNO cycle. In the case of 0.05 v/v_crit the initial production of ¹²C by 3α is delayed, which in turn affects energy production by the CNO cycle. The 0.1 v/v_crit star also encounters the ΩΓ limit at log (L/L_⊙) = 8.256.

Finally, for the 9000 M_⊙ models shown in panel d of Figure 1), the trends are also quite similar to those of the 1000 M_⊙ models but the difference in luminosities in the fastest versus slowest rotating models is even smaller, log (L/L_⊙) = 0.007. At such high masses, all models reach the red side of the HR diagram before the end of core hydrogen burning, with near-zero equatorial velocities. The value of Γ_edd is also quite similar for the three rotation rates, about 0.965. The 0.05 v/v_crit model is on the verge of the ΩΓ limit and 0.1 v/v_crit model encounters it at log (L/L_⊙) = 8.514.

3.2. Rotation and transport of chemical species

We now consider the impact of diffusion coefficients on stellar structure at different stages of evolution. The impact of rotation will be most evident in the 0.1 v/v_crit models, so we show the evolution of these coefficients in the 1000, 5000 and 9000 M_⊙ stars at this rate in Figure 2. All panels a–i depict the inner zones of these models with focus on the boundary of convective core and radiative envelope.

Fig. 2. Evolution of diffusion coefficients in the 0.1 v/vcrit models. The lavender and off-white regions show convective and radiative zones, respectively. Multiple diffusion coefficients overlap a given zone, but the dominant coefficient is the one with the highest magnitude along the y-axis. Dconv is the convective diffusion coefficient, Kther is the thermal diffusivity coefficient, Dshear is the shear (vertical) mixing coefficient, Dh is the horizontal mixing coefficient, and Deff is the effective mixing coefficient. Top row: 1000 M⊙ model at the onset, middle, and end of core hydrogen burning. Centre row: 5000 M⊙ model at the beginning, middle, and end of core hydrogen burning. Bottom row: 9000 M⊙ model at the start, middle, and end of core hydrogen burning.

At the start of core hydrogen burning in panels a, d, g, we find that D_conv is the dominant mixing coefficient in the convective core and the mixing of chemical species in this zone is assumed to be nearly instantaneous. The radii of the convective cores for all three stars shows an intuitive trend, where the 9000 M_⊙ model has the largest (29 R_⊙), followed by 5000 M_⊙ (17.5 R_⊙) and 1000 M_⊙ (9 R_⊙). The thermal diffusivity coefficient, K_ther, is responsible for damping temperature perturbations on a timescale l²/K. In the context of shear instabilities, a large K_ther lowers the effective Richardson number by smoothing out buoyancy forces. Hence the “secular shear” develops when thermal diffusion is fast compared to the turnover of shear eddies. In context of these models, higher K_ther (more efficient radiative diffusion) enhances mixing by weakening stabilizing thermal gradients, and this is evident in radiative zones for all three models throughout the main sequence.

At the boundary of the convective core with a strong mean molecular weight gradient, ∇_μ, we find that the horizontal mixing coefficient D_h is the strongest in magnitude. As seen in Equation (3), D_h is inversely proportional to D_eff, the effective mixing coefficient, but in spite of this dependence it is D_eff that determines chemical transport at the boundary of convective core. This is because D_eff is directly proportional to the radial amplitude of the meridional circulation velocity, and circulation is driven by the thermal imbalance (baroclinicity) created by centrifugal distortion and the resulting horizontal temperature gradients. These temperature gradients are strongest at the convective – radiative boundary, which leads to a much larger meridional circulation velocity than D_h that results in chemical transport dominated by D_eff. Moving outward in the radiative regions, the chemical transport is governed by the vertical mixing coefficient D_shear all the way to the surface of all models. All three models at this evolutionary stage have identical stellar structures and all the mixing coefficients behave in the same manner at all masses.

The stellar structures of the three models begin to show differences at the middle of core hydrogen burning. For the 5000 and 9000 M_⊙ stars (panels e, and h in Figure 2), intermediate convective zones are formed starting at radii of 20 and 30 R_⊙, respectively. These intermediate convective zones are formed as a result of the radiative temperature gradient, ∇_rad, exceeding the adiabatic temperature gradient, ∇_ad, in the local region. The increase in ∇_rad is likely due to the change in local opacity, energy generation rate, and mean molecular weight profile. As in the convective core, chemical transport in these intermediate convective zones is dominated by D_conv. Every other mixing coefficient trend is identical to that at the onset of core hydrogen burning.

At the end of core hydrogen burning (panels c, f, and i in Figure 2), the convective core in all models has receded and helium is being accumulated in the core. The 5000 and 9000 M_⊙ stars are now dominated by intermediate convective zones all the way to the surface, making these models 70–90% convective. At this stage in evolution, transport of chemical species is largely done by D_conv. Conclusively, these results show that D_eff-driven meridional circulation dominates chemical transport at the convective–radiative boundary. Subsequent radiative regions are mixed by D_shear which yields a more uniform envelope composition that can influence burning lifetimes and core growth. Since most of the stellar structure is convective, the transport of chemical species is dominated by the physics of convection.

3.3. Rotation and transport of angular momentum

Rotation also impacts the transport of angular momentum in these stars, but not directly via the diffusion coefficients. In these models (see Eq. (4)), the dominant mechanism for angular-momentum transport is the advective flux associated with meridional circulation, quantified by the radial velocity component U₂(r). Shear instabilities, through the diffusion coefficient D_shear, provide a complementary, diffusive redistribution of angular momentum. Horizontal turbulence, D_h, acts to enforce nearly uniform angular velocity on isobars but contributes minimally to the net transport. Together, these processes establish a balance between advection and diffusion, enabling efficient angular-momentum redistribution in SMSs.

3.4. Supermassive stars and the ΩΓ limit

In order to quantify the effect of the ΩΓ-limit on our SMS models, we focus first on the set with the highest initial rotation (v/v_crit = 0.1), which provides an upper bound on rotationally induced effects. Figure 3 (top panel) shows the evolution of the surface spin parameter Ω/Ω_crit as a function of the local Eddington factor Γ. All three masses start the main sequence at substantial fractions of break-up (Ω/Ω_crit ≃ 0.17, 0.26, 0.33 at Γ ≃ 0.79, 0.93, 0.96 respectively) and rapidly spin up. The value for Ω/Ω_crit reaches ≃0.7 when the ages of the 1000, 5000, and 9000 M_⊙ stars are 2.8 × 10⁴ yr, 7.0 × 10⁵ yr, and 8.2 × 10⁵ yr, respectively, beyond which the centrifugal acceleration g_rot rivals the reduced gravity g_grav(1 − Γ) and mechanical winds must ensue (e.g. Langer 1998). After peaking (up to Ω/Ω_crit ≃ 0.87 and 0.74 for the two highest masses), all models spin down as mass and angular momentum are shed.

Fig. 3. Effects of the ΩΓ limit on the 1000, 5000, and 9000 M⊙ models at 0.1 v/vcrit. Top: Evolution of the surface spin parameter Ω/Ωcrit vs. Eddington factor Γ. The horizontal dashed line marks Γ = 0.639, above which the ΩΓ-limit (Eq. (7)) becomes relevant. Bottom: ratio of radiatively reduced to classical critical velocity, vcrit, 2/vcrit, 1, vs. stellar age. The decline from the initial values of ∼0.69, 0.40, and 0.32 to minima of ∼0.38, 0.25, and 0.19 quantifies the suppression of break-up speed by radiative acceleration. The shaded zones represent the age where the ΩΓ-limit is active.

Figure 3 (bottom panel) plots the ratio v_crit, 2/v_crit, 1 as a function of age. This ratio follows

$$\begin{array}{c}\hfill \frac{v_{\mathrm{crit},2}}{v_{\mathrm{crit},1}}=\frac{3}{2}\sqrt{1-\mathrm{\Gamma }(t)}\frac{R_{\mathrm{eb}}}{R_{\mathrm{e}}},\end{array}$$

and falls from ∼0.69, 0.40, 0.32 at zero-age to minima of ∼0.38, 0.25, 0.19 near each model’s maximum Γ (Eqs. (6)–(7)). The more massive stars, with Γ_max ≳ 0.96, suffer the strongest reduction of their effective break-up speed, demonstrating that even moderate rotation forces them well into the ΩΓ-limit regime and amplifies mass-loss rates by orders of magnitude (Maeder & Meynet 2000). To determine a lower boundary in rotation for the ΩΓ limit, we repeated this analysis for the v/v_crit = 0.05 runs. We find that models in the 8000–10 000 M_⊙ range cross Ω/Ω_crit ≃ 0.7 during the main sequence as well, indicating that mechanical mass loss will set in even at this rotation rate.

In summary, the ΩΓ-limit imposes a stringent cap on SMS rotation: any approach to Ω/Ω_crit ≳ 0.7 triggers strong, mechanical mass loss that prevents sustained high-spin states. This regulatory mechanism, initially proposed for massive stars over a range of metallicities applies equally to supermassive Pop III stars, constraining their angular momenta and influencing their ultimate fates (Maeder & Meynet 2000). However, the extent of mass loss needed to surpass this limit is quite low in numerical simulations. For instance, the 9000 M_⊙ star at 0.1 v/v_crit only loses 147 M_⊙ over its life. This mass, although significant, is only 1.6% of its total mass. To better quantify mass loss, full 3D hydrodynamic simulations for mass-loss outflows from SMSs are required.

3.5. Estimates on the core rotation rates

Figure 4 complements our HR-diagram discussion by showing the mass–weighted rotation velocity of the convective core, v_c, from ZAMS to the end of core H burning for our three Pop III masses. Each panel contains the tracks that start with v/v_crit = 0.01,  0.05, and 0.10; the colour scale gives the instantaneous core radius R_c. For an initial rotation of v/v_crit = 0.1, the 1000, 5000, and 9000 M_⊙ models decelerate almost linearly from initial core rotational velocity ≃160, 210, and 250 km s⁻¹, respectively, to 40–70, 50–80, and 60–90 km s⁻¹ by the end of core H burning. Metal-free winds are negligible, so the loss of core angular momentum must be mediated by internal transport: meridional advection extracts J from the fully convective region and feeds it into the differentially rotating envelope, overcoming the modest contraction of R_c. The decline rate steepens with mass because the moment of inertia grows faster when larger fractions of the star participate in convection.

Fig. 4. Convective-core rotation velocity vs. age for models at three rotation rates. The second y-axis shows the radius of the core in solar radii. Left: 1000 M⊙ models, Centre: 5000 M⊙ models, Right: 9000 M⊙ models.

At fixed mass the vertical spacing of the three tracks reflects the imposed surface spin. After an adjustment phase lasting ∼(0.1–0.15) t_MS the slopes converge, signalling that the core has attained near solid-body rotation and evolves thereafter under the same angular-momentum sink, with a similar rate of change of core velocity with time (dv_c/dt). A global fit to the combined data gives $v_{\mathrm{c}}\propto R_{\mathrm{c}}^{0.3\pm 0.1}$, consistent across all masses and rotation rates and confirming that the decay of J_c dominates over structural changes.

The 5 × 10³ M_⊙, 0.10 v/v_crit model deviates briefly from the other stars from t ≃ 5.5 × 10⁵ and 9 × 10⁵ yr: a short-lived convective shell at r ≃ 25 R_⊙ drains angular momentum from the core, flattens the v_c(t) slope, and disappears once the local μ-gradient is re-established. Model diagnostics show that many intermediate convective shells form at the core-envelope interface; the shell taps angular momentum from the core, flattens the v_c(t) curve, and disappears once the local μ gradient is restored, after which the track rejoins the general trend.

Hydrogen burning occupies ≃90% of the lifetime. Assuming the linear v_c(t) trend extends through the remaining nuclear stages, we combine the late-MS v_c with the measured convective-core masses M_c to estimate the Kerr parameter of the direct-collapse black hole (DCBH). Table 2 gives results for two limits: (i) the DCBH accretes only the envelope above the CO core (M_BH = M_c ) and (ii) the entire star collapses (M_BH = M_⋆ ). The most extreme model, 9000 M_⊙ with a core velocity v_c ≃ 120 km s⁻¹, still reaches only a spin parameter a_*^(core) ≈ 0.05, far below the Thorne limit. Full collapse reduces the spin to a_*^(full) ≈ 0.02. If future simulations yield higher spins, the excess angular momentum must be gained after core-H exhaustion, for example by renewed shear during He burning, secular contraction in the final ∼10⁴ yr, or late-stage fallback.

The spin estimate conserves angular momentum, not linear speed: $J_{\mathrm{c}}=k{M}_{\mathrm{c}}{R}_{\mathrm{c}}^2(v_{\mathrm{c}}/R_{\mathrm{c}})$ is carried unchanged into the collapse, while the physical radius shrinks by five to six orders of magnitude. The Kerr parameter $a_{\ast }=c{J}_{\mathrm{c}}/(G{M}_{\mathrm{BH}}^2)$ therefore remains fixed even though the horizon rotates at v​ = ​a_*c. Any higher IMBH spin must be accumulated after core-H exhaustion, for example through shear generated in core-He burning, secular contraction during the final 10⁴ yr, or late fallback.

4. Discussion

4.1. Comparison with non-rotating models

We compare our rotating SMS tracks, specifically the v/v_crit = 0.1 models, with the non-rotating Pop III models of Nandal et al. (2025a) at the same X_c (Fig. 5). The convective core mass fraction (M_cc/M_tot) reveals how rotation alters core growth against steep radiative and compositional gradients. At 1000 M_⊙, rotation increases the mid-MS core fraction to M_cc/M_tot = 0.780 compared to 0.744 without rotation (+4.6%), and the end-MS core reaches 0.521 versus 0.512 (+1.8%). During the early MS, meridional circulation (U₂) dredges hydrogen across the μ-gradient, sustained by shear instabilities (D_shear), which enlarges the convective core. As the envelope inflates and shear weakens late in the MS, core growth under rotation is partly counteracted, limiting the end-MS gain.

Fig. 5. Mcc/Mtot vs. central hydrogen Xc for non-rotating (dash-dot) and v/vcrit = 0.1 (solid) models at 1000 (blue), 5000 (green), and 9000 (red) M⊙.

At 5000 M_⊙, rotating models show a mid-MS M_cc/M_tot = 0.741 versus 0.738 non-rotating (+0.4%), but by end-MS the core grows from 0.555 to 0.570 (+2.7%). The higher luminosity steepens the radiative gradient (∇_rad), yet rotation amplifies vertical shear mixing (D_shear ∼ 10¹² cm² s⁻¹) and effective diffusion (D_eff ∼ 10⁹ cm² s⁻¹), allowing more hydrogen to penetrate the core throughout the MS. Reduced envelope density late in evolution further enhances this mixing.

At 9000 M_⊙, the rotating mid-MS core fraction (0.736 at X_c ≈ 0.113) falls below the non-rotating 0.774 (−4.9%) because near-Eddington envelopes (Γ_edd ≳ 0.96) initially suppress meridional flows. As the star evolves, centrifugal support reduces the effective gravity (g_eff), meridional velocities recover, and by H exhaustion the rotating core reaches 0.585 versus 0.541 non-rot (+8.1%). This late-MS core boost highlights how rotation can overcome radiative damping to sustain mixing in the most massive SMS.

These end-MS core enhancements translate into MS lifetime extensions of ∼7%, 14%, and 20% at 1000, 5000, and 9000 M_⊙, respectively. At matched X_c, rotating models are cooler by Δlog T_eff ≈ −0.76, −0.04, and ∼0, and luminosities shift by Δlog L ≈ +0.04, −0.03, and +0.002. Eddington factors at these points decrease from 0.894 to 0.864, 0.977 to 0.961, and 0.993 to 0.965, respectively, affecting the onset of mechanical mass loss.

4.2. Comparison with other works

Lee & Yoon (2016) evolved one-dimensional, accreting Pop III protostars at Ṁ ≃ 10⁻¹ M_⊙ yr⁻¹ and showed that once M ≈ 5–7 M_⊙, envelope inflation drives Γ → 1 so that g_eff(1 − Γ) = 0 is reached at v/v_crit ≪ 1, halting further accretion by enforcing the ΩΓ–limit and capping growth at ∼20–40 M_⊙ while keeping radii compact (R ≲ 50 R_⊙). They interpret this as a complete termination of protostellar evolution beyond that mass. Our supermassive models begin at 10³–10⁴ M_⊙ with initial v/v_crit = 0.01–0.10, well below the ∼0.7 threshold for ΩΓ–driven break-up. Although our SMSs cross the ΩΓ–limit, they shed ≲2% of their mass yet continue core H burning to X_c → 0, demonstrating that this combined radiative-centrifugal limit regulates but does not terminate stellar evolution once formation is complete.

Yoon et al. (2012) computed a grid of rotating Pop III stars from 10 to 1000 M_⊙, including magnetic torques and rotational mixing. They find that chemically homogeneous evolution (CHE) occurs for 13 ≤ M ≤ 84 M_⊙ at v/v_crit ≳ 0.4 and ceases above ∼190 M_⊙ when envelope inflation and angular-momentum loss suppress mixing. Although the CHE threshold parallels the ΩΓ–limit onset at moderate masses, our SMS models at v/v_crit ≤ 0.10 remain outside the CHE regime but nonetheless exhibit significant core-mass enhancements via meridional circulation and shear mixing (Sect. 4.1), illustrating that sub-CHE rotation can still measurably alter supermassive stellar structure.

With the formalism developed by Maeder & Meynet (2000), Tsiatsiou et al. (2024) recently used GENEC to compute fast-rotating (v/v_crit = 0.7) Pop III models in the 9–120 M_⊙ range. They found that despite approaching the ΩΓ limit, their models lost ≤1% of their mass. Our SMS models differ in that they are much more massive and luminous, but are rotating much more slowly (v/v_crit ≤ 0.1). The key distinction is that the high Eddington factors in our models (Γ_Edd ≥ 0.96) drastically lower the true critical velocity (v_crit, 2), as seen in the bottom panel of Figure 3. This means that even modest rotation is sufficient to push the star against the ΩΓ limit. Like Tsiatsiou et al. (2024), we find that this results in regulation via moderate mass loss (≤2%) rather than catastrophic disruption, confirming that this physical process operates across a wide range of stellar masses.

Haemmerlé et al. (2018) followed accretion and rotation in Pop III protostars up to supermassive scales and showed that the ΩΓ limit requires the accreted angular momentum to be ≲1% of the Keplerian value. This yields surface velocities v/v_crit ≲ 0.1–0.2 to avoid premature break-up and allow growth to 10⁵ M_⊙. The constraint here enforces slow rotation and negligible centrifugal deformation during accretion. Our rotating SMS models lie precisely within the slow rotation window identified by Haemmerlé et al. (2018) - and although they cross the ΩΓ limit during core H burning, they shed ≲2% of their mass and exhaust hydrogen to X_c → 0, demonstrating that the same ΩΓ regulation applies post-assembly without halting stellar evolution.

Finally, it is worth considering the effect of metallicity. Our models are strictly for metal-free (Pop III) stars. The formation of SMSs may also be possible in very low-metallicity gas (Chon & Omukai 2020). The presence of even trace amounts of metals would significantly enhance mass-loss rates from line-driven winds. As shown in recent work (Sibony et al. 2024), this would alter the angular momentum evolution and could change the final fate and black hole remnant properties, representing an important avenue for future research.

5. Conclusion

We have performed the first simulations of rotating supermassive Pop III star evolution with GENEC. These models examine the interplay of rotation-induced mixing, angular momentum transport, and the ΩΓ limit to advance our understanding of 1000–10 000 M_⊙ primordial stellar evolution. Our main conclusions can be summarized as follows:

• Rotation notably extends the main-sequence lifetime of supermassive Pop III stars by approximately 7% at 1000 M_⊙, 14% at 5000 M_⊙, and 20% at 9000 M_⊙ relative to non-rotating models. This is primarily due to enhanced chemical mixing that enlarges the convective core by up to 8.1% in the most massive cases.

• The ΩΓ limit significantly constrains stellar rotation, triggering mechanical mass loss once surface rotation approaches approximately 70% of the critical velocity. However, total mass loss (not accounting for pulsational mass loss) remains modest (≲2%), suggesting regulation rather than cessation of stellar evolution.

• Internal angular momentum transport via meridional circulation and shear instabilities efficiently redistributes angular momentum, substantially decelerating core rotation. Predicted dimensionless spin parameters for resulting DCBHs remain low (a_* ≤ 0.05 for core collapse scenarios). This implies that rotating SMSs are natural progenitors of slowly rotating black hole seeds. If rapidly spinning seeds are required by observations, they must acquire additional angular momentum through post-main-sequence processes, such as late-stage accretion or fallback.

• Chemical mixing in rotating SMSs is influenced by convective mixing within core and intermediate shells, meridional circulation-driven effective mixing (D_eff ∼ 10⁹ cm² s⁻¹), and shear-induced vertical mixing (D_shear ∼ 10¹² cm² s⁻¹), each critically shaping chemical gradients and core evolution. However, as the structure of the stars is primarily convetive, D_conv is the dominant mixing coefficient for chemical transport throughout the main sequence.

In comparison with previous non-rotating Pop III models, rotation consistently produces larger convective cores and longer stellar lifetimes, a result consistent with lower-mass massive stars. Unlike earlier lower-mass Pop III studies encountering the ΩΓ limit, our supermassive models demonstrate only moderate mass loss and continued stable evolution through core hydrogen exhaustion. These findings align with recent theoretical insights that rotation significantly impacts supermassive star structure without necessarily triggering catastrophic mass loss.

Future investigations will extend these rotating stellar models into higher rotation rates and advanced stages of nuclear burning (Griffiths et al. 2025), explicitly calculating nucleosynthetic yields essential for interpreting high-redshift abundance observations. Additionally, implementing metallicity-dependent mass-loss prescriptions will enhance realism and applicability to early-universe conditions. Inclusion of magnetic instabilities (Spruit 2002; Fuller et al. 2019; Griffiths et al. 2022) will further refine angular momentum transport modelling, which might be particularly important in light of the observations made for low-mass stars for which asteroseismic constraints are available (Aerts et al. 2019; Buldgen & Eggenberger 2023). Finally, exploring the combined effects of rotation and accretion from the pre-main-sequence phase will provide comprehensive evolutionary pathways crucial for understanding the formation and ultimate fate of the earliest supermassive stars.

Acknowledgments

D.N. would like to thank Prof. André Maeder for all the valuable discussions and inspiration for this work. GB acknowledges fundings from the Fonds National de la Recherche Scientifique (FNRS) as a postdoctoral researcher. JR acknowledges support from the Royal Society and Research Ireland through the University Research Fellow programme under grant number URF/R1/191132 and acknowledges support from the Research Ireland Laureate programme under grant number IRCLA/2022/1165. T.E.W. acknowledges the support of the Canadian Space Agency (CSA) [23EXPROSS2] and the support of the Natural Sciences and Engineering Research Council of Canada (NSERC).

References

All Tables

All Figures

Fig. 1. SMS evolution on the HR diagram. The colourbar indicates central H mass fractions, Xc, and the shaded red regions show the isoradii lines: the lightest shade corresponds to log (R/R⊙) = 1 and the darkest shade to log (R/R⊙) = 4. The white outline around the stellar tracks correspond to the evolutionary stage where the ΩΓ-limit is encountered. Top row: 1000–10 000 M⊙ stars at initial rotation rates v/vcrit = 0.01 (a), 0.05 (b), and 0.10 (c). Bottom row: comparison of rotation rates at a given mass, 1000 M⊙ (d), 5000 M⊙ (e), and 9000 M⊙ (f).

In the text

Fig. 2. Evolution of diffusion coefficients in the 0.1 v/vcrit models. The lavender and off-white regions show convective and radiative zones, respectively. Multiple diffusion coefficients overlap a given zone, but the dominant coefficient is the one with the highest magnitude along the y-axis. Dconv is the convective diffusion coefficient, Kther is the thermal diffusivity coefficient, Dshear is the shear (vertical) mixing coefficient, Dh is the horizontal mixing coefficient, and Deff is the effective mixing coefficient. Top row: 1000 M⊙ model at the onset, middle, and end of core hydrogen burning. Centre row: 5000 M⊙ model at the beginning, middle, and end of core hydrogen burning. Bottom row: 9000 M⊙ model at the start, middle, and end of core hydrogen burning.

In the text

Fig. 3. Effects of the ΩΓ limit on the 1000, 5000, and 9000 M⊙ models at 0.1 v/vcrit. Top: Evolution of the surface spin parameter Ω/Ωcrit vs. Eddington factor Γ. The horizontal dashed line marks Γ = 0.639, above which the ΩΓ-limit (Eq. (7)) becomes relevant. Bottom: ratio of radiatively reduced to classical critical velocity, vcrit, 2/vcrit, 1, vs. stellar age. The decline from the initial values of ∼0.69, 0.40, and 0.32 to minima of ∼0.38, 0.25, and 0.19 quantifies the suppression of break-up speed by radiative acceleration. The shaded zones represent the age where the ΩΓ-limit is active.

In the text

	Fig. 4. Convective-core rotation velocity vs. age for models at three rotation rates. The second y-axis shows the radius of the core in solar radii. Left: 1000 M⊙ models, Centre: 5000 M⊙ models, Right: 9000 M⊙ models.
In the text

	Fig. 5. Mcc/Mtot vs. central hydrogen Xc for non-rotating (dash-dot) and v/vcrit = 0.1 (solid) models at 1000 (blue), 5000 (green), and 9000 (red) M⊙.
In the text