© D. Gagnier et al. 2019

1. Introduction

The evolution of the rotation rate of stars is one of the open challenges of current stellar physics. The rotation rate of a star indeed depends on several un-mastered magneto-hydrodynamic mechanisms. The most important ones may be those that transport and/or extract angular momentum within the stellar interior and at the surface of the star, and in the first place the losses due to radiation-driven winds, possibly modified by the presence of a magnetic field. Angular momentum losses depend on various phenomena but in particular on the mass loss distribution at the surface of the star. It is clear that a strong mass loss at the equator of the star is more efficient at extracting angular momentum than a strong mass loss at the pole. Moreover, the shape of a fast rotating star strongly deviates from the spherical symmetry and its spheroidal shape emphasises the anisotropy of the wind.

Until now, stellar evolution codes cope with this question using more or less sophisticated recipes. Often, the spherical symmetry of the models is coupled to a very simplified modelling of mass and angular momentum loss where the star is pealed off at a given rate (e.g. Woosley et al. 1993; Ekström et al. 2012). This rate may depend on some general parameters of the star (luminosity, effective temperature, etc.).

In the present work we make a step forward in the modelling of the rotational evolution of stars by using the 2D ESTER models that have been designed by Espinosa Lara & Rieutord (2013) and Rieutord et al. (2016). These models describe the 2D steady structure of a rapidly rotating star. They can also be used to investigate the rotational evolution of early-type stars along the main sequence, using a simple method that we implement here to compute the change in the hydrogen mass fraction in the convective core of the model, and which provides an acceptable description of the main sequence evolution. With these state-of-the-art 2D models we have now access to the latitude variations of surface quantities that matter for mass loss, namely effective gravity and effective temperature, as the star evolves. We can thus investigate in some details the consequences of the ensuing anisotropic radiation-driven winds and in particular the effects of a mass-flux jump at some effective temperature, as has been suggested by Vink et al. (2001).

In a preliminary study, Gagnier et al. (2019, hereafter referred to as Paper I) have revisited the concept of critical angular velocity for stars with strong surface radiative acceleration. In this study, we first used the ω-model of Espinosa Lara & Rieutord (2011) to derive an analytical expression for the critical angular velocity. Briefly, the ω-model assumes that radiative flux and effective gravity are anti-parallel and that the stellar mass is centrally condensed. We show that up to 90% of the critical angular velocity, the ω-model flux does not diverge more than 10% from the flux given by a full 2D-model. Our new expression of the critical angular velocity differs from that of Maeder & Meynet (2000) pioneer work, and turns out to be very close to the Keplerian angular velocity for all MS evolutionary stages, at least for stellar models of mass less that 40 M_⊙. We explain this small difference by the combined effects of gravity darkening and reduced equatorial opacity caused by the centrifugal force, and conclude that the standard Keplerian angular velocity remains a very good approximation for critical angular velocity. Gagnier et al. (2019) also designed a local mass-flux (ṁ) and angular momentum flux prescription based on the modified CAK theory (Pauldrach et al. 1986) and calibrated with Vink et al. (2001) 1D models. The discontinuity included in our ṁ − T_eff relation (the so-called bi-stability jump) makes the mass-flux of rapidly rotating star controlled by either a single-wind or a two-wind regime (respectively SWR and TWR). In the SWR, we find a maximum mass-flux at the poles, while in the TWR both mass and angular momentum losses are strongly enhanced in equatorial regions. This 2D mass loss prescription now opens the door to the present study of the main-sequence evolution of rotation in rapidly rotating early-type stars.

This paper is organised as follows. In Sect. 2 we present the ESTER code and describe the nuclear time evolution framework. In Sect. 3 we concentrate on the rotational evolution of early-type stars assuming that mass and angular momentum losses are negligible, which is certainly valid for intermediate-mass stars up to ∼7 M_⊙ at metallicity close to solar. In Sect. 4 we consider evolution with mass loss and the associated angular momentum loss, which is appropriate for more massive star evolution. Conclusions follow in Sect. 5.

2. The ESTER code

The ESTER code self-consistently computes the steady state of an isolated rotating star. The models include the 2D axisymmetric structure and the large-scale flows (differential rotation and meridional circulation) driven by the baroclinicity of radiative regions. A short description of the ESTER models can be found in Paper I and more details in Espinosa Lara & Rieutord (2013) and Rieutord et al. (2016).

2.1. A simplified scheme for temporal evolution on the main sequence

As far as time evolution is concerned, we shall concentrate in this study on the main sequence and assume that the star remains in a quasi-steady state. We thus assume that the drivers of evolution are weak enough that all adjustments to the steady state occur on timescales that are much shorter than the ones imposed by the drivers, namely nuclear burning and mass loss.

Since we consider only early-type stars, nuclear reactions are located in the convective core, which is fully mixed. Focusing on the evolution of rotation, we do not need a detailed network of nuclear reactions and simplify the evolution of the mass fraction of hydrogen X by considering its relation to the nuclear energy production ε_*. Heat being essentially produced by the transformation of protons into ⁴He nuclei, we assume that, locally,

$$\begin{array}{c}\hfill \frac{\partial X}{\partial t}+\mathit{v}·\mathbf{\nabla }X=-\frac{4m_{\mathrm{p}}{\epsilon }_{\ast }}{Q},\end{array}$$(1)

where m_p is the proton mass,

$$\begin{array}{c}\hfill Q=(4m_{\mathrm{p}}-m_{\mathrm{He}})c^2\simeq 4.3\times {10}^{-12}J,\end{array}$$(2)

is the energy released by the fusion of the four protons, and ε_* = ε_pp + ε_CNO is the nuclear energy production per unit mass from both pp-chain and CNO cycle.

In order to get an expression for dX_core/dt, we need to calculate the rate of change of the mass of hydrogen in the convective core. As shown in Appendix A, we have

$$\begin{array}{c}\hfill \frac{\mathrm{d}{X}_{\mathrm{core}}}{\mathrm{d}t}=-\frac{4m_{\mathrm{p}}{\overline{\epsilon }}_{\ast }}{Q},\end{array}$$(3)

where ${\overline{\epsilon }}_{\ast }$ is the mass average of ε_*, namely

$$\begin{array}{c}\hfill {\overline{\epsilon }}_{\ast }=\frac{\int {\epsilon }_{\ast }\mathrm{d}M}{M_{\mathrm{core}}}=\frac{{\int }_0^{\pi }{\int }_0^{{\eta }_{\mathrm{core}}}{\epsilon }_{\ast }\rho r_{\zeta }{r}^{2}sin\theta \mathrm{d}\zeta \mathrm{d}\theta }{{\int }_0^{\pi }{\int }_0^{{\eta }_{\mathrm{core}}}\rho r_{\zeta }{r}^{2}sin\theta \mathrm{d}\zeta \mathrm{d}\theta }=\frac{L_{\mathrm{core}}}{M_{\mathrm{core}}}·\end{array}$$(4)

L_core and M_core are respectively the luminosity and the mass of the convective core. η_core is the core fractional radius, (ζ, θ) are the spheroidal coordinates used to describe the star (its surface is at ζ = 1) and r_ζ = ∂r/∂ζ comes from the Jacobian associated with spheroidal coordinates (e.g. Rieutord et al. 2016). Finally, the discretised time evolution of X_core is given by

$$\begin{array}{c}\hfill X_n^{\mathrm{core}}=X_{n-1}^{\mathrm{core}}-\mathrm{\Delta }t\frac{4m_{\mathrm{p}}}{Q}{\overline{\epsilon }}_{\ast },\end{array}$$(5)

where n refers to the time-step of the computations. We now have a “clock” in the code allowing us to monitor the evolution of X_core with time. In addition, we simplify the radial profile of the mean molecular weight at the outer boundary of the convective core by imposing a step-function.

To check that our simplifications provide an acceptable evolution (especially in terms of lifetime on the MS), we compare the time variation of the core hydrogen mass fraction of the simplified ESTER set-up with the more realistic modelling used in the stellar evolution Geneva code. Figure 1 shows the evolution of X_core/X₀ at Z = 0.02 with initial chemical mixture of Grevesse & Noels (1993), as a function of time for a non-rotating 5 M_⊙ star computed with both the ESTER and the Geneva code (Ekström, priv. comm.), where X₀ is the initial hydrogen mass fraction on the zero age main sequence (ZAMS). The models are computed without convective overshoot (this is true for all models in the paper), and without mass loss. The result is that the time at which X_core vanishes, is slightly overestimated by the ESTER code. The difference is typically 20%. We accept such a difference since it is not our purpose to give quantitative predictions, but to rather exhibit the main qualitative features of fast rotation evolution in early-type stars. In the rest of the paper, we use X_core/X₀ as a proxy of time evolution on the MS.

Fig. 1. Evolution of Xcore as a function of time for models of 5 M⊙, Z = 0.02 (GN composition, Grevesse & Noels 1993) computed with ESTER using the simple scheme for hydrogen burning (Eq. (5)) and with the Geneva code (Ekström, priv. comm.). These are non-rotating models computed with no mass loss and no overshooting.

2.2. The quasi-steady state approximation

As mentioned above, the star can be considered in quasi-steady state if its relevant timescales are shorter than the timescales imposed by the drivers. Here, the relevant timescale is the Eddington-Sweet timescale T_ES, which corresponds to the time required for the redistribution of angular momentum (i.e. for baroclinic modes to be damped out, see Rieutord 2006) and defined as

$$\begin{array}{c}\hfill T_{\mathrm{ES}}=T_{\mathrm{KH}}\frac{\mathit{GM}}{{\mathrm{\Omega }}_{\mathrm{eq}}^2{R}^3},\end{array}$$(6)

where Ω_eq is the equatorial angular velocity, and T_KH = GM²/(RL) is the Kelvin-Helmholtz timescale. The comparison between the nuclear and mass loss timescales, T_nucl and T_ml, and the Eddington-Sweet timescale at ZAMS for ESTER 2D-models of 5, 10 and 20 M_⊙ stars initially rotating at ω_i = Ω_eq, i/Ω_k = 0.3 where Ω_k is the Keplerian angular velocity is given in Table 1. It shows that the Eddington-Sweet timescale is at least one order of magnitude shorter than the nuclear timescale and at least two orders of magnitude shorter than the mass loss timescale for our simulations, meaning that the angular momentum has enough time to be redistributed during the stellar evolution, and that the quasi-steady state is a good approximation.

3. Evolution at constant mass and angular momentum

As a first step, we concentrate on “low-mass” early-type stars throughout the Main Sequence (MS). Indeed, for such stars, we can neglect mass and angular momentum losses. This is a good approximation for M ≲ 7 M_⊙ (this limit actually depends on the considered value of metallicity). We thus enforce constant angular momentum and constant mass, namely

$$\begin{array}{c}\hfill L_z={{\displaystyle \int }}_{(V)}{r}^2{sin}^2\theta \mathrm{\Omega }\rho \mathrm{d}V=\mathrm{Cst}\mathrm{and}M={{\displaystyle \int }}_{(V)}\rho \mathrm{d}V=\mathrm{Cst}.\end{array}$$

We compute 2D quasi-steady models with the ESTER code and follow the nuclear evolution according to Eq. (5).

3.1. Evolution of surface stellar parameters

We first consider the main sequence evolution of surface stellar parameters of a 5 M_⊙ model at two different metallicities Z = 0.02 and Z = 10⁻³ (typical of Population II stars) and initially rotating at 50% of the critical angular velocity. Because such stars exhibit very weak radiative acceleration, their critical angular velocity is derived from the Ω-limit, namely when the centrifugal acceleration balances gravity at equator (e.g. Paper I). In that case, the critical angular velocity is the Keplerian angular velocity ${\mathrm{\Omega }}_{\mathrm{k}}=\sqrt{GM/R_{\mathrm{eq}}^3}$. When this critical angular velocity is reached, the spin-up has to stop and the star loses angular momentum such that it remains below and near critical rotation. This angular momentum loss is assumed to occur through mechanical mass transfer into a decretion disc (e.g., Okazaki 2004; Krtička et al. 2011).

To describe the evolution of the rotational properties of the models, we focus on the evolution of the equatorial and polar radii, R_eq and R_p, the surface flattening ϵ = 1 − R_p/R_eq, the ratio of the equatorial angular velocity to the Keplerian one ω, and the linear equatorial velocity V_eq = R_eqΩ_eq. Figures 2–5 display the evolution of these quantities along the MS.

Fig. 2. Equatorial and polar radius as a function of the fractional abundance of hydrogen in the convective core for a 5 M⊙ star initially rotating at ωi = 0.5 for Z = 0.02 and Z = 10−3. The stellar radius in case of no rotation has been added for comparison in full and dashed black lines for Z = 0.02 and 10−3 respectively.

3.1.1. Stellar radius and surface flattening

The mean stellar radius is well known to increase along the MS, the model with the lowest Z being more compact (e.g. Yoon et al. 2006; Ekström et al. 2011). Figure 2 shows that rotating stars expand faster at the equator than at the pole. Indeed, the centrifugal effect affects only weakly the polar regions, and does not contribute to the polar radius growth; however, it does affect the equator significantly. The different behaviour of polar and equatorial radii can be seen in Fig. 3 where the flattening of the star ϵ is plotted along the main sequence evolution. The rapid expansion of the equatorial radius yields a rapid growth of ϵ for the two values of Z. The stellar expansion due to nuclear evolution is slower for the less metallic star, since its keplerian angular velocity Ω_k decreases more slowly, leading to a slower growth of ω and thus of ϵ.

Fig. 3. Same as Fig. 2, but for the surface flattening ϵ = 1 − Rp/Req.

3.1.2. Surface rotation

Figure 4 shows the evolution of the angular velocity ratio ω = Ω_eq/Ω_k throughout the MS. The fact that it increases as X_core/X₀ decreases can be easily explained. Indeed, assuming slow rotation, thus weak flattening and weak differential rotation, we can estimate the moment of inertia of the star as

Fig. 4. Same as Fig. 2, but for the angular velocity at the equator in units of the equatorial critical angular velocity ω = Ωeq/Ωk.

$$\begin{array}{c}\hfill I={\displaystyle \int }{r}^2\mathrm{d}m=AM{R}^2,\end{array}$$(7)

where A informs us on the stellar mass distribution. The smaller A, the more centrally condensed the star. Figure 6 shows the evolution of A = I/MR² as a function of the stellar radius along the MS for a 5 M_⊙ star initially rotating at ω_i = 0.1 and with Z = 0.02. For this model, a rough fit shows that A decreases as ∼R^−0.76. Therefore, as the star evolves along the MS, it becomes more and more centrally condensed. At constant angular momentum L = AMR²Ω = Cst, and as the star expands, Ω decreases as R^−1.24. Ω thus decreases more slowly than Ω_k R^−1.5, leading to an increased angular velocity ratio ω throughout the MS. The foregoing remark shows that the natural trend of the evolution of rotation is an increase of ω, namely an increase of criticality along the MS. The full calculation with 2D models confirms this behaviour at high rotation rate (Fig. 4). Metallicity also plays a role in the evolution of ω: as mentioned before, a less metallic star expands more slowly (see Fig. 2) and thus Ω_k decreases less rapidly, leading to a slower increase of ω. We stress that the star does not rotate faster as it evolves through the MS even though ω increases. Actually, it is the opposite up to the so-called hook at the end of the MS, as shown by Fig. 5. After the hook, nuclear heat generation no longer compensates the radiated energy leading the star to an overall contraction (Iben 1967). The conservation of angular momentum then leads to an increase of equatorial angular velocity Ω_eq.

Fig. 5. Same as Fig. 2, but for the equatorial velocity and angular velocity.

Fig. 6. MS evolution of the ratio A = I/(MR2) as a function of the stellar radius for a 5 M⊙ star initially rotating at ωi = 0.1 and with Z = 0.02. A rough fit represented with a black line gives A R−0.76.

3.1.3. Gravity darkening

The flattening of rapidly rotating stars gives rise to gravity darkening, which reduces the flux at the equator (von Zeipel 1924). Figure 7 shows the variation of the ratio between the local effective temperature and the polar effective temperature as a function of colatitude for the 5 M_⊙ model with Z = 0.02, at ZAMS and for various values of ω_i. The latitudinal variation of this ratio at Z = 10⁻³ is not shown because almost identical to the case with Z = 0.02, with a relative difference never exceeding 0.3%. The effective temperature drop between pole and equator is thus almost exactly the same for the two metallicities, as expected from Fig. 3 showing that both models have the same surface flattening at ZAMS.

Fig. 7. Variation of the effective temperature at ZAMS as a function of colatitude for a 5 M⊙ star at Z = 0.02 for various values of ωi.

3.2. Initial angular velocity requirements for criticality

To appreciate the initial states (at ZAMS) that lead to the critical velocity before the end of the MS, we compute a grid of models with constant masses 5 ≤ M/M_⊙ ≤ 10 and with Z = 0.02 and Z = 10⁻³, backward in time. Hence we start with models rotating near criticality (here we take ω = 0.9, i.e. quasi-criticality) at various hydrogen mass fractions in the core (X_core, c/X₀), or equivalently at various stage along the MS, and increase the hydrogen mass fraction until X_core = X₀, namely until ZAMS is reached.

Figure 8 shows the ZAMS angular velocity ratio ω_i needed to reach quasi-criticality during the MS for various hydrogen mass fractions in the core at quasi-criticality, as a function of the stellar mass M/M_⊙, when Z = 0.02. The more massive the star, the smaller ω_i is needed to reach critical rotation on the MS, in particular for X_core/X₀ ≤ 0.5. Moreover, the smaller the fractional abundance of hydrogen left in the convective core at quasi-criticality, the smaller ω_i needs to be. Indeed, as mentioned before, stellar expansion tends to make the angular velocity ratio ω increase as the star evolves through the MS. Therefore, the smaller the fractional abundance of hydrogen left in the convective core at quasi-criticality, the longer the evolution along the MS and thus the more ω has increased. At ZAMS, the star therefore does not need to be rotating very rapidly to reach criticality near the end of the MS.

Fig. 8. Values of the ZAMS angular velocity ratio ωi required to reach quasi-criticality (ω = 0.9) at different values of the hydrogen mass fraction in the core, as a function of the stellar mass M/M⊙ for Z = 0.02, and without any mass loss. These predictions result from computations backward in time.

We now take a look at the effect of lower metallicity on the initial conditions required to reach criticality before the end of the MS. The comparison between the ZAMS angular velocities at equator in units of the equatorial critical angular velocity as a function of the mass for Z = 0.02 and Z = 10⁻³ is shown in Fig. 9. We clearly see that for a lower metallicity, the value of the initial angular velocity required to reach criticality at a given hydrogen content on the MS is higher than at higher metallicity. This reflects the slower expansion of less metallic stars.

Fig. 9. Same as Fig. 8 for Xcore, c/X0 = 0.9 and 0.7, and for Z = 0.02 and 10−3.

4. Evolution with angular momentum loss

We now consider stars more massive than ∼7 M_⊙ for which mass loss from radiation-driven winds can no longer be neglected.

4.1. The critical angular velocity at the ΩΓ-limit

The critical angular velocity of massive stars is expected to differ from the Keplerian angular velocity because of radiative acceleration (see Maeder & Meynet 2000). Using the θ-dependent radiative flux expression from the ω-model of Espinosa Lara & Rieutord (2011), we have shown in Paper I that, in this approximation, the angular velocity becomes critical when the rotation-dependent Eddington parameter at equator

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{\Omega }}(\pi /2)={\mathrm{\Gamma }}_{\mathrm{eq}}{(1-\frac{{\mathrm{\Omega }}^2{R}_{\mathrm{eq}}^3}{\mathit{GM}})}^{-2/3},\end{array}$$(8)

reaches unity, where

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{eq}}=\frac{\kappa (\pi /2)L}{4\pi cGM},\end{array}$$(9)

is the standard Eddington parameter at equator. The critical angular velocity then reads

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{c}}={\mathrm{\Omega }}_{\mathrm{k}}\sqrt{1-{\mathrm{\Gamma }}_{\mathrm{eq}}^{3/2}}.\end{array}$$(10)

Thus, it is expected that the critical angular velocity is reduced by radiative acceleration at the stellar surface. However, as shown in Paper I, this reduction is actually quite small because of the combined effects of gravity darkening and opacity reduction at equator. We further illustrate this property in Fig. 10 where we plot the evolution of Γ_eq and Γ_Ω(π/2) for a 15 M_⊙ ESTER 2D-model. Unlike the no-mass loss case discussed in Paper I, mass loss prevents the star from reaching criticality (see below) making Γ_Ω(π/2) always finite. For such a star, Γ_eq ≲ 0.27, implying that Ω_c/Ω_k ≳ 0.93 all the time. This result shows that the difference between critical and Keplerian angular velocity remains small. Thus, we shall continue to express Ω_eq as a fraction of the equatorial Keplerian angular velocity Ω_k to measure the distance to criticality.

Fig. 10. Γeq and ΓΩ(π/2) as a function of the fractional abundance of hydrogen in the convective core for a 15 M⊙ ESTER 2D-model initially rotating at ωi = 0.5, with and without mass loss (respectively solid and dashed lines). The black line shows the evolution of Γeq for ω = 0. The two vertical lines delineate the two phases of the two-wind regime (see Sect. 4.5).

4.2. Mass and angular momentum loss

Similarly to the case of evolution without angular momentum loss, we evolve the star with ESTER through the MS by decreasing the fractional abundance of hydrogen in the convective core X_core/X₀. However, at each step, we also remove some mass and angular momentum off the star, for that we use the clock given by the time evolution of X_core in Eq. (3).

We estimate the mass-loss rate using the local mass-flux prescription of Paper I, which can be seen as a local equivalent of the modified CAK theory (Pauldrach et al. 1986) for rotating winds. Although not clearly assessed by observations (see the discussion in Paper I), we also include the bi-stability jump in our local mass loss modelling.

The total mass and angular momentum loss rates for a rotating star read

$$\begin{array}{c}\hfill \dot{M}={\displaystyle \int }\dot{m}(\theta )\mathrm{d}S\mathrm{and}\dot{\mathcal{L}}={\displaystyle \int }\dot{\ell }(\theta )\mathrm{d}S.\end{array}$$(11)

The area element at the stellar surface is

$$\begin{array}{c}\hfill \mathrm{d}S=2\pi R^2(\theta )\sqrt{1+\frac{R_{\theta }^2}{R^2(\theta )}}sin\theta \mathrm{d}\theta ,\end{array}$$

where R(θ) is the θ-dependent radius of the star and R_θ = ∂R/∂θ (Rieutord et al. 2016). From Paper I, the local mass-flux is estimated via

$$\begin{array}{cc}\hfill \dot{m}(\theta )& =\frac{4}{9}\frac{\alpha (\theta )k\prime {(\theta )}^{1/\alpha \prime (\theta )}}{v_{\mathrm{th}}(\theta )c}\hfill \\ \hfill & \times {\left[\frac{c}{{\kappa }_{\mathrm{e}}(1-\alpha )}\left(|g_{\mathrm{eff}}(\theta )|-\frac{{\kappa }_{\mathrm{e}}F(\theta )}{c}\right)\right]}^{\frac{\alpha \prime (\theta )-1}{\alpha \prime (\theta )}}\hfill \\ \hfill & \times F{(\theta )}^{1/\alpha \prime (\theta )},\hfill \end{array}$$(12)

and the local angular momentum flux reads

$$\begin{array}{c}\hfill \dot{\ell }(\theta )=\dot{m}(\theta )\mathrm{\Omega }(\theta )R({\theta )}^2{sin}^2\theta ,\end{array}$$(13)

where Ω(θ) is the θ-dependent surface angular velocity. In Eq. (12), g_eff(θ) is the local gravito-centrifugal acceleration, F(θ) is the local radiative flux, κ_e is the mass absorption coefficient for electron scattering, and v_th ≡ (2k_BT_eff/m_Fe)^1/2 is the local iron thermal velocity. α and k are the CAK force multiplier parameters (FMPs), and α′=α − δ where δ is another FMP related to the effect of radial changes of ionisation when moving outward in the wind. We assume δ not to vary on the stellar surface, and equal to 0.1, a typical value for hot stars at metallicity close to solar (Abbott 1982). We also use the parametrisation of α and k at Z = 0.02 from Paper I, namely

$$\begin{array}{c}\hfill \alpha (T_{\mathrm{eff}})=\{\begin{array}{cc}0.45\hfill & \mathrm{if}{T}_{\mathrm{eff}}\le 10\mathrm{kK}\hfill \\ 1.5\times {10}^{-5}{T}_{\mathrm{eff}}+0.3\hfill & \mathrm{if}10\mathrm{kK}<T_{\mathrm{eff}}\le 20\mathrm{kK}\hfill \\ 5\times {10}^{-6}{T}_{\mathrm{eff}}+0.5\hfill & \mathrm{if}20\mathrm{kK}<T_{\mathrm{eff}}\le 40\mathrm{kK}\hfill \\ 0.7\hfill & \mathrm{if}{T}_{\mathrm{eff}}>40\mathrm{kK},\hfill \end{array}\end{array}$$(14)

and

$$\begin{array}{c}\hfill \begin{array}{ccc}& k(T_{\mathrm{eff}})\simeq \hfill & \hfill \{\begin{array}{cc}exp(-2.15\times {10}^{-4}{T}_{\mathrm{eff}}+2.41)\hfill & \mathrm{if}{T}_{\mathrm{eff}}\le 20\mathrm{kK}\hfill \\ -3.00\times {10}^{-6}{T}_{\mathrm{eff}}+0.22\hfill & \mathrm{if}20\mathrm{kK}<T_{\mathrm{eff}}\le T_{\mathrm{eff}}^{\mathrm{jump}}\hfill \\ 1.16\times {10}^{-6}{T}_{\mathrm{eff}}+0.08\hfill & \mathrm{if}{T}_{\mathrm{eff}}>T_{\mathrm{eff}}^{\mathrm{jump}}.\hfill \end{array}\end{array}\end{array}$$(15)

$T_{\mathrm{eff}}^{\mathrm{jump}}$ is the effective temperature of the bi-stability jump, that is, the temperature below which Fe IV recombines into Fe III, and is defined as (Vink et al. 2001)

$$\begin{array}{c}\hfill T_{\mathrm{eff}}^{\mathrm{jump}}=61.2+2.59⟨\rho ⟩,\end{array}$$(16)

where ⟨ρ⟩ is the characteristic wind density at 50% of the terminal velocity of the wind. It is given by

$$\begin{array}{c}\hfill log⟨\rho ⟩=-14.94+3.2{\mathrm{\Gamma }}_{\mathrm{e}},\end{array}$$(17)

where

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{e}}=\frac{{\kappa }_{\mathrm{e}}L}{4\pi cGM},\end{array}$$(18)

is the classical Eddington parameter. In 1D, the Fe IV–Fe III recombination theoretically leads to an increased global mass-loss rate and a decreased terminal velocity of the wind. The effective temperature at the surface of rotating stars is however θ-dependent and the bi-stability jump can therefore occur locally if there is a latitude on the stellar surface where $T_{\mathrm{eff}}=T_{\mathrm{eff}}^{\mathrm{jump}}$. We note that Eq. (15) has been obtained assuming that the global mass-loss rate in the non-rotating regime Ṁ = 4πR²ṁ is equivalent to the one of Vink et al. (2001) (some caveats on the use of this mass loss prescription for the calibration of k are discussed in Paper I). k has been calibrated with ESTER models at ZAMS, although it slightly varies with the evolutionary stage along the MS due to the evolution of g_eff. We find that Δk = k(X_core/X₀ = 0.1) − k(X_core/X₀ = 1)> 0 never exceeds 20% of k(X_core/X₀ = 1) for all T_eff. Since our minimum value for α′ is 0.35 and ṁ k^1/α′, the local mass-flux may be underestimated by a factor ∼1.7, at most. This difference is acceptable considering the other approximations used and the uncertainties of the wind model.

4.3. A remark on mass loss effects without nuclear evolution

At this stage, it is useful to recall how the radius of a star changes when a small amount of mass is removed. The outcome is well-known for the case of fully convective stars, which are well represented by a n = 1.5 polytrope. For such stars, if a small amount of mass is removed from the surface, the remaining mass expands (Chandrasekhar 1967). On the other hand, radiative stars, which are well represented by n >  3 polytropes, contract when subject to mass loss (Heisler & Alcock 1986). This is also true for ESTER 2D-models of massive stars with a radiative envelope, as can be seen in Fig. 11 showing the evolution of the polar and equatorial radius of a 15 M_⊙ star initially rotating at ω_i = 0.5, from which we have extracted a total mass of ΔM ≃ 0.4 M_⊙, keeping X_core/X₀ = 1 constant.

Fig. 11. Equatorial and polar radius as a function of the mass for a 15 M⊙ star with Z = 0.02 initially rotating with ωi = 0.5 and losing a total amount of mass of ΔM ≃ 0.4 M⊙.

Unfortunately, the interpretation of the ω-dependence on mass and angular momentum loss is rather complicated due to the non-linearity of the processes involved. While the angular momentum loss of massive stars only affects the evolution of ω by strengthening the decrease of Ω_eq, the effect of slower expansion due to mass loss is twofold, since it leads to a slower decrease of both Ω_k and Ω_eq. Whether mass loss itself tends to accelerate or slow down the natural trend of ω to increase throughout the MS is therefore surely model dependent. Figure 12 shows the evolution of both Ω_eq and Ω_k with and without mass and angular momentum loss, for a 15 M_⊙ ESTER 2D-model with Z = 0.02 and for ω_i = 0.5. For this specific model, we see that the evolution of ω almost only depends on the evolution of Ω_k, that is because the effect of mass loss and angular momentum loss on the evolution of Ω_eq roughly compensate.

Fig. 12. Angular velocity (red) and equatorial Keplerian angular velocity (black) as a function of the fractional abundance of hydrogen in the convective core Xcore/X0 for a 15 M⊙ star with a ZAMS angular velocity ratio ωi = 0.5, with and without mass loss (respectively solid and dashed lines).

4.4. The two different wind regimes

As shown in Paper I, the radiation-driven wind of rotating massive stars can be in two distinct regimes. The first one is a SWR, characterised by a maximum mass-flux at the poles, which is effective when the effective temperature of the star is either greater or smaller than $T_{\mathrm{eff}}^{\mathrm{jump}}$ for all θ. The second one is a TWR, which appears when there exists a colatitude θ_jump on the stellar surface where $T_{\mathrm{eff}}({\theta }_{\mathrm{jump}})=T_{\mathrm{eff}}^{\mathrm{jump}}$. The latter regime is characterised by a maximum mass-flux at equator and a stronger global mass and angular momentum loss rates. Intermediate-mass stars (or well-evolved massive stars, i.e. late B- and A-supergiants) can be in a cold-sided SWR, that is, with $T_{\mathrm{eff}}<T_{\mathrm{eff}}^{\mathrm{jump}}$ for all θ, characterised by stronger global mass and angular momentum loss rates than in a TWR. Finding the general rules that govern the evolution of rotating massive stars turned out to be quite cumbersome since numerous particular cases pop up. Hence, to get an idea of the rotational evolution of massive stars, we focus on the case of a typical 15 M_⊙ star initially rotating at 50% of the Keplerian angular velocity.

4.5. The example of a 15 M_⊙ 2D-model with ω_i = 0.5

Figures 13 and 14 respectively show the MS evolution of the mass-loss rate, the angular velocity ratio, and the equatorial effective temperature for an ESTER 2D-model with M = 15 M_⊙, Z = 0.02, and ω_i = 0.5. We find the star to start its MS evolution in a (hot-sided) SWR, that is, with an equatorial effective temperature larger than $T_{\mathrm{eff}}^{\mathrm{jump}}$, until X_core/X₀ ≃ 0.4. When X_core/X₀ ≲ 0.4, T_eff, eq is smaller than the effective temperature of the bi-stability jump and the star enters a TWR, then remains in this regime for the rest of the MS. In the SWR, log Ṁ increases roughly linearly and ω increases similarly to the case with no mass loss, which is described in Sect. 3. The TWR can be divided into two phases of evolution. A first phase for 0.4 ≲ X_core/X₀ ≲ 0.3 (Phase 1), where Ṁ rapidly increases while ω keeps increasing but more slowly. And a second phase for X_core/X₀ ≲ 0.3, where log Ṁ goes back to a roughly linear increase, while ω decreases. In Fig. 15 we illustrate the MS evolution of the global mass-loss rate against the corresponding surface-averaged effective temperature ${\overline{T}}_{\mathrm{eff}}$ of the model. As expected (see Paper I), we find the TWR to be reached even if the surface-averaged effective temperature is larger than $T_{\mathrm{eff}}^{\mathrm{jump}}$. It underlines that the average T_eff is not an appropriate quantity to determine the wind regime. Snapshots of the different phases along the evolution for this model are shown in Fig. 16. We now give some explanation for the evolution of rotation and mass-loss rate in the three phases that we have identified.

Fig. 13. Top: evolution of the mass-loss rate Ṁ (in M⊙ yr−1) as a function of the fractional abundance of hydrogen in the convective core Xcore/X0 for a 15 M⊙ star with Z = 0.02 and for ωi = 0.5. The vertical dashed line marks the evolutionary time at which the star reaches a TWR. The vertical dotted line marks the transition between Phase 1 and 2. Bottom: same but for the angular velocity at equator in units of the equatorial critical angular velocity ω = Ωeq/Ωk. The red dashed line corresponds to the evolution of ω with constant mass and angular momentum throughout the MS.

Fig. 14. As Fig. 13, but for the evolution of the equatorial effective temperature Teff, eq (in K) as a function of the fractional abundance of hydrogen in the convective core Xcore/X0 for a 15 M⊙ star with Z = 0.02 and for ωi = 0.5. The black full line corresponds to the evolution of the effective temperature of the bi-stability jump $T_{\mathrm{eff}}^{\mathrm{jump}}$.

Fig. 15. Evolution of the mass-loss rate Ṁ as a function of the mean effective temperature ${\overline{T}}_{\mathrm{eff}}$ for a 15 M⊙ ESTER 2D-model with Z = 0.02 and for ωi = 0.5. The star evolves from left to right with a mean effective temperature approaching 30 kK at ZAMS, and ${\overline{T}}_{\mathrm{eff}}\simeq 22$ kK at TAMS. The black curve corresponds to case of a non-rotating 15 M⊙ star with Z = 0.02 for which the mass-loss rate has been calculated using Vink et al. (2001) prescription. The vertical dashed line marks the evolutionary time at which the star reaches a TWR. The vertical dotted line marks the transition between Phase 1 and 2.

Fig. 16. Snapshots of the local mass-flux ṁ in units of the polar mass-flux at ZAMS ṁp,ZAMS (left) and of the local angular momentum flux $\dot{\ell }$ in units of the maximum angular momentum flux at ZAMS ${\dot{\ell }}_{\max ,\mathrm{ZAMS}}$ (right) for a 15 M⊙ ESTER 2D-model with Z = 0.02 and for ωi = 0.5, at four different time steps: Xc ≃ 1, 0.3, 0.03 and 0. The first one corresponds to ZAMS, the second corresponds to the end of Phase 1, the third is the end of Phase 2, and the last snapshot corresponds to the very end of the MS evolution. The star is viewed with an inclination i = 70°.

4.5.1. The single-wind regime

The MS evolution of massive stars is well known to exhibit an increase of stellar luminosity L, and a decrease of surface averaged effective temperature ${\overline{T}}_{\mathrm{eff}}$. This is a consequence of the increase of central temperature and density as hydrogen burning proceeds. From Eqs. (12), (14) and (15), it is thus clear that Ṁ should increase throughout the MS, for all rotation rates, in agreement with mass loss scaling relations (e.g. de Jager et al. 1988; Kudritzki et al. 1995). Moreover, as shown in Paper I, the effect of rotation on Ṁ is relatively weak in the SWR (see also Maeder & Meynet 2000; Petrenz & Puls 2000). The effect of the growth of ω on the evolution of Ṁ is negligible compared to the effects of secular evolution.

Figure 13 (bottom) shows that, similarly to the case with no mass loss, ω increases, but more slowly. That is to be expected since angular momentum losses tend to emphasise the decrease of Ω_eq as evolution proceeds. Additionally, mass loss also slows down the stellar expansion as mass loss itself tends to make the star contract (see Sect. 4.3). Ω_k therefore decreases more slowly when mass loss is accounted for (see Fig. 12).

4.5.2. The first phase of the two-wind regime

For our 15 M_⊙ model, the first phase of the TWR phase corresponds to the 0.4 ≲ X_core/X₀ ≲ 0.3 time period. During this phase, the evolution of ω and Ṁ are tightly coupled. Indeed, at X_core/X₀ ≃ 0.4, $T_{\mathrm{eff},\mathrm{eq}}\simeq T_{\mathrm{eff}}^{\mathrm{jump}}$, and the region where the mass-flux is increased due to the local reach of bi-stability limit, is confined to the equator. At this stage, the mass and angular momentum losses are not strong enough to prevent ω from increasing. The growth is only reduced. The two phenomena controlling the evolution of Ṁ are the nuclear evolution of the star and the evolution of ω. In this phase, both tend to make θ_jump migrate poleward. Figure 17 shows the evolution of the surface fraction of the star, ΔS/S, where the effective temperature is lower than $T_{\mathrm{eff}}^{\mathrm{jump}}$, that is, the portion of enhanced mass-flux due to bi-stability. In this phase, ΔS/S increases, thus leading to an increase of the mass and angular momentum loss rates. This leads to a faster decrease of Ω_eq and a slower decrease of the Keplerian angular velocity Ω_k (Fig. 12) compared to models without mass loss. This process lasts until the mass and angular momentum loss rates are sufficient for ω to start decreasing. At this stage, the star enters the second phase of the TWR.

4.5.3. The second phase of the two-wind regime

For this stellar model, the second phase of the TWR phase corresponds to X_core/X₀ ≲ 0.3, that is, from the point where ω starts decreasing. As ω decreases we would expect θ_jump to migrate back to the equator. However, ΔS/S and thus Ṁ still increase during this phase. Obviously, nuclear (secular) evolution effects still dominate over the consequences of the ω decrease. Hence, θ_jump still migrates poleward. As ω decreases in this phase, the surface flattening ϵ also decreases. The effect of secular nuclear evolution becomes predominant over the effect of rotation on the evolution of equatorial opacity. This leads to a new increase of Γ_eq almost until the end of the MS (see Fig. 10).

Fig. 17. Evolution of the surface fraction of the star where the effective temperature is lower than $T_{\mathrm{eff}}^{\mathrm{jump}}$, for a 15 M⊙ star with Z = 0.02 and for ωi = 0.5. The black dashed line marks the evolutionary time at which the star reaches a TWR. The black dotted line marks the threshold between Phase 1 and 2.

Fig. 18. Angular velocity at the equator in units of the equatorial critical angular velocity ω = Ωeq/Ωk as a function of the fractional abundance of hydrogen in the convective core Xcore/X0 for different initial stellar masses and for Z = 0.02, assuming a ZAMS angular velocity ratio ωi = 0.3 (left) and ωi = 0.7 (right).

After the hook near end of the MS, the star contracts leading to an increase of the mean effective temperature, thus to the migration of θ_jump towards the equator. It hence leads to a decrease of Ṁ and to an ultimate increase of ω just before the end of the MS (see Fig. 13).

4.5.4. Other masses

ESTER predictions for a given initial angular velocity ratio ω_i = 0.3 (left) and ω_i = 0.7 (right), and for various stellar masses, are shown in Fig. 18. It turns out that the evolution of ω is a rather complicated function of the mass and the initial angular velocity, at least for the mass range of our models, namely 7 ≲ M/M_⊙ ≲ 15. This complexity comes from the fact that the local effective temperature of such stars is often close to the effective temperature of the bi-stability jump. Hence, the stars are prone to either start their MS evolution in a hot-sided SWR (i.e. $T_{\mathrm{eff}}>T_{\mathrm{eff}}^{\mathrm{jump}}$ for all θ) and possibly reach a TWR followed by a cold-sided SWR (i.e. $T_{\mathrm{eff}}<T_{\mathrm{eff}}^{\mathrm{jump}}$ for all θ), or start in a TWR and reach a cold-sided SWR. Hence, the evolution of ω may drastically vary when initial conditions are changed.

The modelling of more massive stars is certainly more predictable. Indeed, such stars are always in a hot-sided SWR at ZAMS, and above a certain mass, the expansion and evolution of rotation is not sufficient to reach a TWR. Additionally, more massive stars lose a lot more mass and angular momentum because of their stronger radiation-driven wind, and thus may exhibit a decreasing ω all along the MS. Unfortunately, for numerical reasons, the ESTER code presently does not allow us to study the nuclear evolution of models with masses much higher than ∼20 M_⊙.

5. Conclusions

This work presents a study of the rotational evolution of rapidly rotating early-type stars with the 2D code ESTER that we have updated to follow the decrease of core hydrogen content as a proxy for time evolution along the MS. We first investigated the rotational evolution of intermediate-mass stars with masses in-between 5 and 7 M_⊙. The evolution of such stars along the MS is relatively simple because at first order it can be modelled with constant angular momentum, that is, without considering any mass loss. We have shown that, because of stellar expansion as well as the redistribution of mass in the stellar interior, critical rotation can be reached before the end of the MS provided the initial angular velocity ratio ω_i = Ω_eq, i/Ω_k is high enough. The minimum value of ω_i required to reach critical rotation is smaller for the more massive stars. In this mass range, no star reaches the critical angular velocity during the MS if its initial angular velocity is less than 50% critical.

At masses M ≥ 7 M_⊙, our models show that the rotational evolution is more complex. Indeed, such stars are subject to radiation-driven winds, inducing mass and angular momentum losses, which substantially impact the evolution of the angular velocity ratio ω. Unfortunately, the magnitude of mass-loss rates obtained with CAK-based prescriptions are rather uncertain, because of loosely constrained parameters and the neglect of wind inhomogeneities (see Paper I for more details). More accurate mass-loss rates could be derived from recent hydrodynamically consistent non-local thermodynamical equilibrium stellar atmosphere models (Sander et al. 2017), but such models will certainly remain 1D for some time. Nevertheless, we used the mass-flux prescription detailed in Paper I, based on the modified CAK theory (Pauldrach et al. 1986), where we include, locally, the bi-stability limit. For now, the existence of this jump is still challenged by the observations (e.g. Markova & Puls 2008), but their conclusions are not clear cut. From the point of view of the models, the bi-stability jump has important consequences and hence cannot be overlooked. Indeed, in 2D-models, the bi-stability limit imposes the existence of two-wind regimes that have rather different effects on the rotational evolution of massive stars. Stars are in a SWR when, at all colatitudes, their effective temperature is either larger or smaller than the effective temperature of the bi-stability jump. This wind regime exhibits a maximum mass-flux at the poles. However, when there is a colatitude θ_jump where the effective temperature equals the effective temperature of the bi-stability jump, the star is in a TWR. In this regime, the mass-flux is maximum in the equatorial region and holds a discontinuity at θ_jump.

The main result of this study is that the rotational evolution of massive stars with masses in-between 7 and 15 M_⊙ is very dependent on their initial conditions at ZAMS, namely their mass and initial angular velocity. Indeed, their effective temperature is rather close to the effective temperature of the bi-stability jump. The initial angular velocity therefore controls whether stellar models in this mass range start their MS evolution in a SWR or in a TWR. The different wind regimes and transitions between these regimes as evolution proceeds, can have very different effects on the rotational evolution of such stars, hence the importance of the ZAMS initial conditions. The most common transition from one regime to another is the transition from a hot-sided SWR to a TWR. We have detailed the rotational evolution of a 2D-ESTER model with M = 15 M_⊙ initially rotating with an angular velocity ratio ω_i = 0.5. The SWR-TWR transition is followed by two phases of evolution of the angular velocity ratio. These two phases correspond to a rapid (and continuous) increase of the mass loss and angular-momentum-loss rates leading to a slower growth of ω (compared to the SWR), followed by the decrease of ω accompanied by a slow increase of the mass and angular momentum loss rates.

Another important result is that, using 2D models, and despite the use of Vink et al. (2001) to calibrate our mass loss prescription in the non-rotating case, we do not find a discontinuity of the global mass-loss rate at the bi-stability limit as in 1D models. Instead, the discontinuity lies on the local mass-flux at the equator, which then migrates poleward as evolution proceeds leading to a continuous increase of the global mass-loss rate in the TWR. Because of the local nature of the bi-stability limit in 2D models, the latter could be reached at various rotation rates and stellar masses, thus various surface-averaged effective temperatures. In this work, we have chosen the effective temperature of the bi-stability limit to be that of Vink et al. (2001), namely roughly 25 kK. Taking a lower value, for instance ∼20 kK from Petrov et al. (2016) would shift the stellar models that may reach a TWR during the MS towards less massive stars. Models that remain in the requisite mass range would reach a TWR later on the MS. The effects of a reduced effective temperature of the bi-stability limit might become particularly important for later evolutionary phases of more massive stars (B-supergiants).

All in all, we find that radiation-driven mass and angular momentum losses are responsible for either a slower growth or a decrease of the angular velocity ratio ω along the MS. This often prevents massive stars from reaching critical rotation before the end of the MS.

The present work calls for new investigations to confirm the present results, in particular in the modelling of the mass losses. This is an old problem, which we here face through the modelling of the force multipliers parameters (see references in Paper I). It is known that these parameters depend on the metallicity of the atmosphere sourcing the wind, but the changes to be applied for modelling a population III star, for instance, are still quite uncertain (e.g. Meynet et al. 2008). The attempt of Georgy et al. (2013) finds that a change in metallicity does not bring any significant change in the rotational evolution of a 15 M_⊙ model, but their models are one-dimensional. We also note that there is strong evidence that massive stars undergo periods of super-Eddington conditions (Quataert et al. 2016) in which they could exhibit continuum-driven winds, sufficiently strong to sustain significant mass loss, even for non-metallic stars (e.g., Owocki et al. 2004, 2017; Smith & Owocki 2006).

Beside these open questions related to the influence of metallicity, the present work also calls for an investigation of the coupling between mass losses and interior flows. Indeed, the rotational evolution of the star implies the presence of meridional flows that carry chemicals from the core to the surface. Presently, the 2D evolution given by ESTER models is a succession of steady states, with steady flows, that are hardly usable to predict the transport of chemicals inside the star. New studies, that will include viscous stresses everywhere in the star, and not only in Ekman layers as presently, are needed to give a full 2D-view of the transport process in a radiative envelope and shed new light on the long standing question of rotational mixing.

Acknowledgments

We are particularly grateful to Joachim Puls for his detailed reading, comments and suggestions on the original manuscript. We thank Georges Meynet and Fabrice Martins for enlightening discussions. We are grateful to Sylvia Ekström for providing information to validate our scheme for main sequence temporal evolution. We thank CALMIP – the computing centre of Toulouse University (Grant 2017-P0107). M. Rieutord acknowledges the strong support of the French Agence Nationale de la Recherche (ANR), under grant ESRR (ANR-16-CE31-0007-01), and of the International Space Science Institute (ISSI) for its support to the project “Towards a new generation of massive star models” lead by Cyril Georgy. F. Espinosa Lara acknowledges the financial support of the Spanish MINECO under project ESP2017-88436-R.

References

Appendix A: Rate of change of the mass fraction of hydrogen in the stellar convective core

In this appendix we calculate the rate of change of the mass fraction of hydrogen in the core. We write the local conservation of hydrogen (which is equivalent to Eq. (1)) as

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}(\rho X)+\mathrm{div}(\rho X\mathit{v})=-\frac{4m_{\mathrm{p}}}{Q}\rho {ϵ}_{\ast }.\end{array}$$(A.1)

In order to get an expression for dX_core/dt, we calculate the rate of change of the mass of hydrogen in the core, namely

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}t}(M_{\mathrm{core}}{X}_{\mathrm{core}})=\frac{\mathrm{d}}{\mathrm{d}t}{{\displaystyle \int }}_{\mathrm{core}}\rho X\mathrm{d}V·\end{array}$$(A.2)

This equation can be re-written using Leibniz’s rule as

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}t}(M_{\mathrm{core}}{X}_{\mathrm{core}})={{\displaystyle \int }}_{\mathrm{core}}\frac{\partial }{\partial t}(\rho X)\mathrm{d}V+{{\displaystyle \int }}_{\partial \mathrm{c}ore}\rho X{\mathit{v}}_{\mathrm{c}}·\mathit{dS},\end{array}$$(A.3)

where v_c is the speed of the core boundary. Using Eq. (A.1) and the fact that we consider the core to be fully mixed, that is, X is homogeneous in the core,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{\mathrm{d}}{\mathrm{d}t}(M_{\mathrm{core}}{X}_{\mathrm{core}})& ={\displaystyle -}\frac{4m_{\mathrm{p}}}{Q}{L}_{\mathrm{core}}\hfill \\ \hfill & +X_{\mathrm{core}}(-{{\displaystyle \int }}_{\mathrm{core}}\mathrm{div}(\rho \mathit{v})\mathrm{d}V+{\int }_{\partial \mathrm{c}ore}\rho {\mathit{v}}_{\mathrm{c}}·\mathit{dS}).\hfill \end{array}\end{array}$$(A.4)

Using Leibniz’s rule again, Eq. (A.4) becomes

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{\mathrm{d}}{\mathrm{d}t}(M_{\mathrm{core}}{X}_{\mathrm{core}})& =-\frac{4m_{\mathrm{p}}}{Q}{L}_{\mathrm{core}}+X_{\mathrm{core}}\frac{\mathrm{d}}{\mathrm{d}t}{{\displaystyle \int }}_{\mathrm{core}}\rho \mathrm{d}V\hfill \\ \hfill & =-\frac{4m_{\mathrm{p}}}{Q}{L}_{\mathrm{core}}+X_{\mathrm{core}}\frac{\mathrm{d}{M}_{\mathrm{core}}}{\mathrm{d}t},\hfill \end{array}\end{array}$$(A.5)

which finally leads to

$$\begin{array}{c}\hfill \frac{\mathrm{d}{X}_{\mathrm{core}}}{\mathrm{d}t}=-\frac{4m_{\mathrm{p}}}{Q}\frac{L_{\mathrm{core}}}{M_{\mathrm{core}}}·\end{array}$$(A.6)

All Tables

All Figures

	Fig. 1. Evolution of Xcore as a function of time for models of 5 M⊙, Z = 0.02 (GN composition, Grevesse & Noels 1993) computed with ESTER using the simple scheme for hydrogen burning (Eq. (5)) and with the Geneva code (Ekström, priv. comm.). These are non-rotating models computed with no mass loss and no overshooting.
In the text

	Fig. 2. Equatorial and polar radius as a function of the fractional abundance of hydrogen in the convective core for a 5 M⊙ star initially rotating at ωi = 0.5 for Z = 0.02 and Z = 10−3. The stellar radius in case of no rotation has been added for comparison in full and dashed black lines for Z = 0.02 and 10−3 respectively.
In the text

	Fig. 3. Same as Fig. 2, but for the surface flattening ϵ = 1 − Rp/Req.
In the text

	Fig. 4. Same as Fig. 2, but for the angular velocity at the equator in units of the equatorial critical angular velocity ω = Ωeq/Ωk.
In the text

	Fig. 5. Same as Fig. 2, but for the equatorial velocity and angular velocity.
In the text

	Fig. 6. MS evolution of the ratio A = I/(MR2) as a function of the stellar radius for a 5 M⊙ star initially rotating at ωi = 0.1 and with Z = 0.02. A rough fit represented with a black line gives A R−0.76.
In the text

	Fig. 7. Variation of the effective temperature at ZAMS as a function of colatitude for a 5 M⊙ star at Z = 0.02 for various values of ωi.
In the text

	Fig. 8. Values of the ZAMS angular velocity ratio ωi required to reach quasi-criticality (ω = 0.9) at different values of the hydrogen mass fraction in the core, as a function of the stellar mass M/M⊙ for Z = 0.02, and without any mass loss. These predictions result from computations backward in time.
In the text

	Fig. 9. Same as Fig. 8 for Xcore, c/X0 = 0.9 and 0.7, and for Z = 0.02 and 10−3.
In the text

	Fig. 10. Γeq and ΓΩ(π/2) as a function of the fractional abundance of hydrogen in the convective core for a 15 M⊙ ESTER 2D-model initially rotating at ωi = 0.5, with and without mass loss (respectively solid and dashed lines). The black line shows the evolution of Γeq for ω = 0. The two vertical lines delineate the two phases of the two-wind regime (see Sect. 4.5).
In the text

	Fig. 11. Equatorial and polar radius as a function of the mass for a 15 M⊙ star with Z = 0.02 initially rotating with ωi = 0.5 and losing a total amount of mass of ΔM ≃ 0.4 M⊙.
In the text

	Fig. 12. Angular velocity (red) and equatorial Keplerian angular velocity (black) as a function of the fractional abundance of hydrogen in the convective core Xcore/X0 for a 15 M⊙ star with a ZAMS angular velocity ratio ωi = 0.5, with and without mass loss (respectively solid and dashed lines).
In the text

Fig. 13. Top: evolution of the mass-loss rate Ṁ (in M⊙ yr−1) as a function of the fractional abundance of hydrogen in the convective core Xcore/X0 for a 15 M⊙ star with Z = 0.02 and for ωi = 0.5. The vertical dashed line marks the evolutionary time at which the star reaches a TWR. The vertical dotted line marks the transition between Phase 1 and 2. Bottom: same but for the angular velocity at equator in units of the equatorial critical angular velocity ω = Ωeq/Ωk. The red dashed line corresponds to the evolution of ω with constant mass and angular momentum throughout the MS.

In the text

	Fig. 14. As Fig. 13, but for the evolution of the equatorial effective temperature Teff, eq (in K) as a function of the fractional abundance of hydrogen in the convective core Xcore/X0 for a 15 M⊙ star with Z = 0.02 and for ωi = 0.5. The black full line corresponds to the evolution of the effective temperature of the bi-stability jump $T_{\mathrm{eff}}^{\mathrm{jump}}$.
In the text

Fig. 15. Evolution of the mass-loss rate Ṁ as a function of the mean effective temperature ${\overline{T}}_{\mathrm{eff}}$ for a 15 M⊙ ESTER 2D-model with Z = 0.02 and for ωi = 0.5. The star evolves from left to right with a mean effective temperature approaching 30 kK at ZAMS, and ${\overline{T}}_{\mathrm{eff}}\simeq 22$ kK at TAMS. The black curve corresponds to case of a non-rotating 15 M⊙ star with Z = 0.02 for which the mass-loss rate has been calculated using Vink et al. (2001) prescription. The vertical dashed line marks the evolutionary time at which the star reaches a TWR. The vertical dotted line marks the transition between Phase 1 and 2.

In the text

Fig. 16. Snapshots of the local mass-flux ṁ in units of the polar mass-flux at ZAMS ṁp,ZAMS (left) and of the local angular momentum flux $\dot{\ell }$ in units of the maximum angular momentum flux at ZAMS ${\dot{\ell }}_{\max ,\mathrm{ZAMS}}$ (right) for a 15 M⊙ ESTER 2D-model with Z = 0.02 and for ωi = 0.5, at four different time steps: Xc ≃ 1, 0.3, 0.03 and 0. The first one corresponds to ZAMS, the second corresponds to the end of Phase 1, the third is the end of Phase 2, and the last snapshot corresponds to the very end of the MS evolution. The star is viewed with an inclination i = 70°.

In the text

	Fig. 17. Evolution of the surface fraction of the star where the effective temperature is lower than $T_{\mathrm{eff}}^{\mathrm{jump}}$, for a 15 M⊙ star with Z = 0.02 and for ωi = 0.5. The black dashed line marks the evolutionary time at which the star reaches a TWR. The black dotted line marks the threshold between Phase 1 and 2.
In the text

	Fig. 18. Angular velocity at the equator in units of the equatorial critical angular velocity ω = Ωeq/Ωk as a function of the fractional abundance of hydrogen in the convective core Xcore/X0 for different initial stellar masses and for Z = 0.02, assuming a ZAMS angular velocity ratio ωi = 0.3 (left) and ωi = 0.7 (right).
In the text