© D. Gagnier et al. 2019

1. Introduction

Among the numerous problems that need to be overcome when stars are modelled, those related to rotation are of particular nature in the frame of classical 1D models because rotation breaks the imposed spherical symmetry. Rotating stars are indeed not only distorted by the centrifugal acceleration, but are also pervaded by large-scale flows that carry chemical elements and angular momentum. The importance of these effects has been appreciated for quite some time now (e.g. Maeder & Meynet 2000, and references therein), and specific modelling simplifications are usually included in 1D stellar evolution codes to reproduce the expected effects of global rotation. For instance, the transport of angular momentum that results from small-scale turbulence and large-scale circulation induced by rotation in radiative zones is inserted in 1D evolution models either as an advection-diffusion process following Zahn (1992), Meynet & Maeder (1997), and Maeder & Zahn (1998) (e.g. Geneva code, Eggenberger et al. 2008; STAREVOL, Decressin et al. 2009, Amard et al. 2016; FRANEC, Chieffi & Limongi 2013; CESTAM, Marques et al. 2013) or as a purely diffusive process (e.g. Kepler, Heger et al. 2000; STERN, Yoon & Langer 2005; MESA, Paxton et al. 2013). The associated transport of chemicals (so-called rotation-induced mixing) is always treated as a diffusive process (as justified by Chaboyer & Zahn 1992).

This modelling of rotation effects is only justified for slow rotators (Zahn 1992). Early-type stars are, however, often considered to be fast rotators. The hypotheses and approximations of current prescriptions are therefore no longer valid for such stars. Be-type stars, for instance, are well known to be fast rotators close to the break-up limit (e.g. Porter & Rivinius 2003; Bastian et al. 2017), that is, the centrifugal force nearly balances gravity at equator. These stars show evidence of mass loss that is associated with their near break-up rotation (Carciofi et al. 2008; Krtička et al. 2011; Rivinius et al. 2013; Georgy et al. 2013; Granada et al. 2013; Granada & Haemmerlé 2014). Furthermore, early-type stars may also be very luminous and therefore have high radiation pressure at their surface. The induced radiation-driven wind is responsible for a significant loss of mass and angular momentum, which notably influences the evolutionary paths of massive stars (e.g. Langer 1998; Vink et al. 2010). Because of gravity darkening (e.g. Espinosa Lara & Rieutord 2011), mass loss from rotating massive stars is expected to be anisotropic (e.g. Owocki et al. 1996; Owocki & Gayley 1997; Pelupessy et al. 2000; Maeder & Meynet 2000; Georgy et al. 2011). These anisotropies in turn affect the evolution of rotation and are likely to play a significant role in the interior dynamics of massive stars (Zahn 1992; Maeder 1999; Lignieres et al. 2000; Lau et al. 2011; Rieutord & Beth 2014).

The treatment of fast rotation thus requires developments beyond the current model approximations, although this approach has been extremely useful to make significant progress in the field (e.g. Maeder & Meynet 2015, and references therein). In this context, the achievement of the first self-consistent 2D models of rapidly rotating early-type stars, worked out by Espinosa Lara and Rieutord (e.g. Espinosa Lara & Rieutord 2013; Rieutord et al. 2016), opens the door to exploring the evolution of fast stellar rotators. Such models are expected to provide new constraints on the internal rotation-induced mechanisms as well as on the radiative and mechanical mass loss (e.g. Krtička et al. 2011), which all significantly affect the different predictions and outputs of the 1D stellar evolution models (e.g. Meynet & Maeder 2000; Maeder & Meynet 2000, 2010; Smith 2014; Meynet et al. 2015; Renzo et al. 2017).

The present work aims at investigating two questions that are prerequisite to the study of the evolution of massive rapidly rotating stars: (i) What is the critical angular velocity of a star when radiative acceleration is significant in its atmosphere? (ii) How do mass and angular momentum loss in massive stars depend on rotation?

This paper is organised as follows. In Sect. 2 we reconsider the question of the critical angular velocity in light of ESTER 2D models and the simplified ω-model of Espinosa Lara & Rieutord (2011). We then revisit the prescription of mass loss in fast-rotating stars. To this end, we first focus on deriving a local mass-flux prescription based on the 1D CAK (Castor et al. 1975) and mCAK (Pauldrach et al. 1986) theories for non-rotating stars (Sect. 3). Next, we generalise this prescription to rotating stars. We compute the latitudinal variations of mass and angular momentum fluxes with ESTER 2D models and discuss the effects of rotation on global losses of mass and angular momentum (Sect. 4). We finally summarise our answers to the questions that motivated this work (Sect. 5).

2. Critical angular velocity and the ΩΓ-limit

2.1. The ΩΓ-limit question

2.1.1. Some context

In stars more massive than ∼7 M_⊙ that are close to solar metallicity, radiative acceleration plays a significant role in the (assumed) hydrostatic equilibrium. Total gravity is usually introduced,

$$\begin{array}{c}\hfill {\mathit{g}}_{\mathrm{tot}}={\mathit{g}}_{\mathrm{eff}}+{\mathit{g}}_{\mathrm{rad}},\end{array}$$(1)

where effective gravity (or gravito–centrifugal acceleration) g_eff is supplemented by radiative acceleration,

$$\begin{array}{c}\hfill {\mathit{g}}_{\mathrm{rad}}=\frac{\kappa \mathit{F}}{c},\end{array}$$(2)

where κ is the flux-weighted opacity¹, which we approximate with the total Rosseland mean opacity, F the radiative flux, and c the speed of light.

The so-called ΩΓ-limit introduced by Maeder & Meynet (2000) is reached when g_tot = 0 somewhere on the stellar surface. It is associated with an actual critical angular velocity Ω_c that is different from the Keplerian angular velocity

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{k}}=\sqrt[]{\frac{\mathit{GM}}{R_{\mathrm{eq}}^3}},\end{array}$$(3)

which is the break-up limit (or the Ω-limit) when radiative acceleration can be neglected in total gravity. In this equation, G is the gravitation constant, M is the stellar mass, and R_eq is its equatorial radius.

The expression of the correct critical angular velocity when the effects of radiation cannot be neglected has been debated lively. For instance, Langer (1997, 1998) suggested that stars close to the Eddington limit have a lower critical angular velocity, while Glatzel (1998) stressed that the Eddington parameter, namely

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

where L is the stellar luminosity, has no effect on the critical rotation because of gravity darkening. In an attempt to clarify the debate, Maeder (1999) and Maeder & Meynet (2000; hereafter referred to as MMM) re-investigated the question and found two roots to the equation g_tot = 0. The first gives the Keplerian angular velocity as the critical angular velocity for Eddington parameters smaller than 0.639. The second root yields a critical angular velocity lower than Ω_k that tends to zero when the rotation-dependent Eddington parameter (see Maeder 1999) tends to unity for Eddington parameters larger than 0.639.

Maeder (1999) based his derivation on the model of von Zeipel (1924), which states that the radiative flux F at some colatitude θ on the surface of a rotating star is proportional to the local effective gravity g_eff. For barotropic stars, this leads to

$$\begin{array}{c}\hfill \mathit{F}=-\rho \chi \frac{\mathrm{d}T}{\mathrm{d}P}{\mathit{g}}_{\mathrm{eff}},\end{array}$$(5)

where χ = 4acT³/(3κρ) is the radiative conductivity. Additionally, assuming solid-body rotation, Zahn (1992) obtained

$$\begin{array}{c}\hfill \rho \chi \frac{\mathrm{d}T}{\mathrm{d}P}=\frac{L(P)}{4\pi G{M}_{\star }},\end{array}$$(6)

where

$$\begin{array}{c}\hfill M_{\star }=M(1-\frac{{\mathrm{\Omega }}^2}{2\pi G{\rho }_{\mathrm{m}}}),\end{array}$$(7)

and ρ_m is the mean density of the star. L(P) is the luminosity outflowing across the isobar of pressure P. Maeder (1999) then wrote the radiative flux in the barotropic case as

$$\begin{array}{c}\hfill \mathit{F}=-\frac{L(P)}{4\pi G{M}_{\star }}{\mathit{g}}_{\mathrm{eff}}.\end{array}$$(8)

In the case of shellular rotation, Ω ≃ Ω(r), and following the work of Zahn (1992), Maeder (1999) linearly developed all quantities around their average on an isobar and found the radiative flux in the baroclinic case

$$\begin{array}{c}\hfill \mathit{F}=-\frac{L(P)}{4\pi G{M}_{\star }}(1-\zeta (\theta )){\mathit{g}}_{\mathrm{eff}}.\end{array}$$(9)

We note that Eq. (8) assumes solid-body rotation, while Eq. (9) corresponds to the case of slow rotation (Zahn 1992). Maeder (1999) noted that ζ(θ)∼10⁻² so that according to this model, the ratio F/g_eff depends only mildly on colatitude.

2.1.2. Interferometric observations

Recent progresses on rotating stars, both observational and theoretical, does not confirm this mild dependence, however. Interferometric observations of several rapidly rotating stars (e.g. Monnier et al. 2007; Zhao et al. 2009; Che et al. 2011; Domiciano de Souza et al. 2014) show that if gravity darkening is modelled by a power law such as

$$\begin{array}{c}\hfill F(\theta )\sim g_{\mathrm{eff}}{(\theta )}^{4\beta },\end{array}$$(10)

then β <  1/4 for all the observed stars. Furthermore, the observed exponents decrease as the rotation rate of the stars increases (e.g. Domiciano de Souza et al. 2014). These results are in line with the predictions of the ESTER 2D models, which match the observations well (Espinosa Lara & Rieutord 2011). ESTER 2D models indeed predict that the relation between flux and effective gravity is not a power law, but can be approximated as such as a first step. Models also show that β ≃ 0.25 at low rotation rates, but that β decreases to 0.13 when rotation approaches criticality. This behaviour has implications on the ΩΓ-limit introduced by Maeder & Meynet (2000). These limitations prompt us to revisit this limit with ESTER 2D models.

2.2. The ω-model

Before using full ESTER 2D models, it is worth considering the problem in light of the simplified ω-model of Espinosa Lara & Rieutord (2011). The general purpose of the ω-model is to describe the latitudinal variations in radiative flux of rotating stars in a simpler way than with a full 2D model. To this end, it is assumed that the flux within the radiative envelope of an early-type star can be written as

$$\begin{array}{c}\hfill \mathit{F}=-f(r,\theta ){\mathit{g}}_{\mathrm{eff}},\end{array}$$(11)

where f(r, θ) is some function of the position to be determined. In this assumption, F and g_eff are assumed to be anti-parallel. Espinosa Lara & Rieutord (2011) showed that this is a rather good approximation because full ESTER 2D models show that the angle between the two vectors never exceeds half a degree, even for the most distorted stars.

In a radiative stellar envelope, the function f(r, θ) can be determined from flux conservation equation, namely ∇ ⋅ F = 0, along with the assumption that the stellar mass is rather centrally condensed, so that the Roche model can be used. This implies that the first-order equation of the flux can be completed by the boundary condition

$$\begin{array}{c}\hfill \underset{r\to 0}{lim}f(r,\theta )=\frac{L}{4\pi GM}·\end{array}$$(12)

The equation for f(r, θ) can then be solved analytically (see Espinosa Lara & Rieutord 2011; Rieutord 2016, for details), with the following result:

$$\begin{array}{c}\hfill f(r,\theta )=\frac{L}{4\pi GM}\frac{{tan}^2\psi (r,\theta )}{{tan}^2\theta },\end{array}$$(13)

where ψ(r, θ) is obtained by solving

$$\begin{array}{c}\hfill cos\psi +lntan(\psi /2)=\frac{1}{3}{\omega }^2{r}^3{cos}^3\theta +cos\theta +lntan(\theta /2).\end{array}$$(14)

In this equation r has been scaled with the equatorial radius R_eq and

$$\begin{array}{c}\hfill \omega =\frac{\mathrm{\Omega }}{{\mathrm{\Omega }}_{\mathrm{k}}}\end{array}$$(15)

where Ω is the angular velocity of the star, which is assumed to be uniform (the case of surface differential rotation has been investigated by Zorec et al. 2017).

At the equator, an analytic expression of f can be obtained,

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

so that the equatorial radiative flux reads

$$\begin{array}{c}\hfill \mathit{F}(R_{\mathrm{eq}},\pi /2)=-\frac{L}{4\pi GM}{(1-{\omega }^2)}^{-2/3}{\mathit{g}}_{\mathrm{eff}}.\end{array}$$(17)

In the slow rotation limit, Eq. (17) can be written

$$\begin{array}{c}\hfill \mathit{F}(R_{\mathrm{eq}},\pi /2)\simeq -\frac{L}{4\pi GM(1-\frac{2}{3}{\omega }^2)}{\mathit{g}}_{\mathrm{eff}}.\end{array}$$(18)

In this limit, where R_eq ≈ R_p, this is identical to the MMM expression, which we now obtain for negligible stellar distortions.

Equations (17) and (18) show an important difference: the exact (within the ω-model) expression (17) shows that the ratio F/g_eff diverges when the Ω-limit ω = 1 is approached, while in its slow rotation approximation (18), F/g_eff remains finite. This is in line with von Zeipel’s law, which is valid at low rotation rates and which states that flux and effective gravity are proportional. This important difference now calls for a new investigation of the ΩΓ-limit.

2.3. Critical angular velocity at the ΩΓ-limit: Ideas from the ω-model

2.3.1. Preliminaries

With the expression of the radiative flux from the ω-model, we can derive the critical angular velocity Ω_c corresponding to the ΩΓ-limit. When this limit is reached, then

$$\begin{array}{c}\hfill {\mathit{g}}_{\mathrm{tot}}={\mathit{g}}_{\mathrm{eff}}+{\mathit{g}}_{\mathrm{rad}}=\mathbf{0}\end{array}$$(19)

somewhere at the surface of the star. As in MMM, we introduce a limiting flux from Eq. (2) and the ΩΓ-limit condition g_tot = 0, namely

$$\begin{array}{c}\hfill {\mathit{F}}_{\lim }=-\frac{c}{\kappa }{\mathit{g}}_{\mathrm{eff}}.\end{array}$$(20)

From this expression, we define the rotation-dependent Eddington parameter Γ_Ω(θ) as the ratio of the actual flux F(θ) obtained with the ω-model, and the limiting flux, namely

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{\Omega }}(\theta )=\frac{F(\theta )}{F_{\lim }(\theta )}=\frac{\kappa (\theta )}{c}f(r=1,\theta ).\end{array}$$(21)

Using Eq. (21) we can rewrite Eq. (1) as

$$\begin{array}{c}\hfill {\mathit{g}}_{\mathrm{tot}}={\mathit{g}}_{\mathrm{eff}}[1-{\mathrm{\Gamma }}_{\mathrm{\Omega }}(\theta )].\end{array}$$(22)

The critical angular velocity Ω_c is reached if somewhere on the stellar surface g_tot = 0, that is, if there is a colatitude where either Γ_Ω(θ)=1 or g_eff(θ)=0.

2.3.2. Critical latitude: the equator

In all 2D models both effective temperature and effective gravity are minimum at the equator (e.g. Espinosa Lara & Rieutord 2013). The solution g_eff(θ) = 0 is therefore always reached first at the equator. We now focus on the Γ_Ω(θ) = 1 solution.

We first observe that Γ_Ω(θ) ∝ κ(θ)f(r = 1, θ) should be an increasing function of co-latitude, at least for (very) rapidly rotating stars. f(r = 1, θ) indeed always increases with θ and diverges at the equator when the Keplerian angular velocity is approached. As mentioned before, in all 2D models the effective temperature is minimum at equator, but it is not straightforward how to predict whether the opacity κ increases or decreases with decreasing T_eff (see Appendix A for an attempt). Still, we expect the opacity to vary on the stellar surface, but much less than f(r, θ) near Keplerian angular velocity. According to the ω-model, and because of the equatorial singularity, it is therefore very likely that the solution Γ_Ω = 1 is always first reached at the equator.

Maeder & Meynet (2000) came to the same conclusion regarding the location of the Γ_Ω(θ) = 1 solution on the surface. However, they traced it back to an opacity effect, assuming that the latter increases with decreasing T_eff and thus is highest at the equator. According to ESTER 2D models, this is not the case for rotating stars (see below).

2.3.3. Unique critical angular velocity

Equation (22) shows that there are two solutions for g_tot = 0, and thus two possible critical angular velocities. However, as we show now, the ω-model removes the g_eff = 0 root for the ΩΓ-limit and thus points to a single critical angular velocity. According to the ω-model at equator, Eq. (1) reads

$$\begin{array}{c}\hfill g_{\mathrm{tot}}(\pi /2)=g_{\mathrm{eff}}(\pi /2)+g_{\mathrm{rad}}(\pi /2),\end{array}$$(23)

where

$$\begin{array}{c}\hfill g_{\mathrm{rad}}(\pi /2)=-\frac{\kappa (\pi /2)L}{4\pi cGM}{(1-{\omega }^2)}^{-2/3}{g}_{\mathrm{eff}}(\pi /2),\end{array}$$(24)

and

$$\begin{array}{c}\hfill g_{\mathrm{eff}}(\pi /2)=-R_{\mathrm{eq}}{\mathrm{\Omega }}_{\mathrm{k}}^2(1-{\omega }^2).\end{array}$$(25)

Here g_eff and g_rad are the radial components of the accelerations (thus positive when outwards). We can then write the equatorial total gravity scaled with ${\mathrm{\Omega }}_{\text{k}}^2{R}_{\text{eq}}$ as

$$\begin{array}{c}\hfill {\stackrel{\sim }{g}}_{\mathrm{tot}}(\pi /2)={\omega }^2+{\mathrm{\Gamma }}_{\mathrm{eq}}{(1-{\omega }^2)}^{1/3}-1,\end{array}$$(26)

where Γ_eq is the standard Eddington parameter evaluated at the equator. From Eq. (24), we see that at the equator, the ratio g_rad/g_eff increases as (1 − ω²)^−2/3 with increasing ω, which also implies that if g_eff approaches 0 when ω → 1, g_rad will also tend to 0 but more slowly. Figure 1 shows the scaled total gravity, effective gravity, and radiative acceleration at the equator as a function of ω with Γ_eq = 0.5. The total gravity at the equator has two zeros; the first root corresponds to Γ_Ω(π/2) = 1, and the second root gives g_eff(π/2) = 0.

Fig. 1. Scaled total gravity, effective gravity, and radiative acceleration at the equator as a function of ω with Γeq = 0.5. The total gravity at the equator has two zeros; the first root corresponds to ΓΩ(π/2)=1, and the second to geff(π/2)=0.

For sub-critical rotation (i.e. g_tot <  0 or |g_eff|> g_rad), the star is gravitationally bound. When we increase ω, the equatorial effective gravity |g_eff| decreases faster than g_rad, to the point where |g_eff| = g_rad (equivalently, Γ_Ω(π/2) = 1), at this point, Ω = Ω_c and g_tot = 0. Increasing ω even more would result in a radiative acceleration that surpasses the effective gravity at equator. When this happens, g_tot >  0 and the star is no longer gravitationally bound up to ω = 1 where the second root is reached. The solution Γ_Ω(π/2) = 1 is therefore always reached before ω = 1, when evolution (say) drives the growth of ω.

2.3.4. Critical rotation given by the ω-model

The main difference between the MMM model and ours, in addition to our unique critical angular velocity, comes from the latitudinal variation of Γ_Ω. In MMM models the latitudinal variations in Γ_Ω come from the latitudinal variations in opacity when we discard the small correcting function ζ(θ). As a consequence, if the surface opacity were constant (e.g. with Thomson opacity of electrons), Γ_Ω would reach unity at all latitudes at the same time when ω increases! The ω-model predicts that the ratio between effective gravity and radiative flux depends on latitude and diverges at the equator when criticality approaches. Even in the extreme case of a constant surface opacity, only a small equatorial region therefore becomes unbound at criticality. The ω-model shows that opacity variations over the stellar surface are unimportant for determining the latitude where g_tot = 0 because of the equatorial singularity. This discussion demonstrates that the use of a constant ratio between the surface flux and the effective gravity (the von Zeipel law) as done in the MMM model has an important consequence for determining a critical rotation because it removes the equatorial singularity of the ratio $T_{\text{eff}}^4/g_{\text{eff}}$.

In line with the ω-model and the maximum of Γ_Ω at equator, the condition giving the critical angular velocity Ω_c is

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{\Omega }}(\pi /2)=\frac{\kappa (\pi /2)L}{4\pi cGM}{(1-\frac{{\mathrm{\Omega }}_{\mathrm{c}}^2}{{\mathrm{\Omega }}_{\mathrm{k}}^2})}^{-2/3}=1,\end{array}$$(27)

or equivalently,

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

These equations show that Ω_c is reduced with increasing Eddington parameter compared to Ω_k. Maeder & Meynet (2000) came to the same conclusion, but with a different expression for critical angular velocity, namely ${\mathrm{\Omega }}_{\mathrm{c}}\propto {\mathrm{\Omega }}_{\mathrm{k}}\sqrt{1-{\mathrm{\Gamma }}_{\mathrm{eq}}}$. Because Γ_eq ≤ 1, their ratio critical to Keplerian angular velocity is lower for the same equatorial Eddington parameter. When radiative acceleration effects are weak at the equator, that is, when Γ_eq ≪ 1, we find Ω_c ≃ Ω_k, as expected. This is also the solution of Maeder & Meynet (2000) for critical angular velocity in this regime.

2.3.5. Some conclusions from the ω-model

The analysis based on the ω-model underlines three important points:

1. Formally, the Keplerian angular velocity is never reached. The critical angular velocity such that the centrifugal acceleration overcomes the sum of the gravitational and radiative accelerations at some place on the stellar surface is always lower than the Keplerian angular velocity.

2. This balance of forces is always first reached at the equator when Ω/Ω_k increases because the ratio $T_{\text{eff}}^4/g_{\text{eff}}$ at the equator diverges when criticality is reached.

3. The use of the von Zeipel law, which assumes the proportionality of the surface flux with the effective gravity gives a critical latitude that depends on the latitudinal variations of the surface opacity and is therefore not necessarily located at the equator.

These conclusions based on the ω-model immediately raise the question of the accuracy of this model. This is the next point that we discuss in light of observations and full 2D ESTER models.

2.4. The ΩΓ-limit with ESTER 2D models

2.4.1. Interferometric observations (again)

We first briefly return to observations. As described above (Sect. 2.1.2), interferometric observations of fast-rotating stars all show that β <  1/4 when the surface flux distribution is assumed to vary as $T_{\text{eff}}\propto g_{\text{eff}}^{\beta }$. Moreover, they clearly show that β decreases with increasing rotation (Domiciano de Souza et al. 2014). From this discussion we can now interpret the fact that 4β <  1 and decreases with increasing rotation as evidence for a divergence of the ratio $T_{\text{eff}}^4/g_{\text{eff}}$ at the equator when criticality is approached.

2.4.2. Accuracy of the ω-model

The assumptions of the ω-model are that the flux vector is anti-parallel to the effective gravity, the gravitational field is that of a point mass (the Roche model), and the rotation is uniform. The last two of these approximations probably entail the largest errors. They can be appreciated by comparing the flux latitudinal distribution of the ω-model with the output of 2D ESTER models. A comparison has been made in Espinosa Lara & Rieutord (2011), but here we focus on the relative difference between the flux of the two models.

We computed the flux from the ω-model, Eq. (11), where L, M, ω, r, and g_eff were taken from the output of ESTER 2D models. For two 2D ESTER zero-age main sequence (ZAMS) models of 15 M_⊙ and 40 M_⊙, we computed the relative difference between the fluxes of the ESTER and ω-model, namely

$$\begin{array}{c}\hfill \frac{\delta F}{F}=\frac{|F_{\mathrm{ESTER}}-F_{\omega }|}{F_{\mathrm{ESTER}}}\end{array}$$(29)

as a function of co-latitude. The result is shown in Fig. 2.

Fig. 2. Relative difference between the radiative flux of the ω-model and the ESTER model as a function of co-latitude for a 15 M⊙ (top) and a 40 M⊙ (bottom) ZAMS-star with various angular velocity ratios. The black line corresponds to an evolved 15 M⊙ ESTER 2D model with ω = 0.9 and a fractional abundance of hydrogen in the convective core Xcore/X0 = 0.5. X0 is the initial hydrogen mass fraction at ZAMS. The minimum of each curve corresponds to a sign change of FESTER − Fω.

On the ZAMS, and for rotation rates of up to 90% of the Keplerian angular velocity, the relative difference between the fluxes remains lower than 10%. For angular velocity ratios lower than 50%, this difference drops to less than one percent, making the ω-model quite reliable for most of the rapidly rotating stars.

This comparison has been made at ZAMS. As stars evolve along the MS, they become more and more centrally condensed (this is discussed in the follow-up paper) and thus better satisfy the Roche approximation. Therefore, the relative deviation between the radiative flux of ESTER 2D models and the analytic ω-model probably never exceeds 10% for stars with ω ≤ 0.9 during the MS. This is illustrated in Fig. 2 (top) with a 15 M_⊙ model at mid-MS rotating with ω = 0.9. Clearly, the relative difference is reduced compared to the ZAMS model.

2.4.3. The ΩΓ-limit with ESTER 2D models

The current ESTER 2D models describe the steady state of a rotating star with a convective core and a radiative envelope, that is, an early-type star. Compared to previous attempts of making stellar models in two dimensions (e.g. Roxburgh 2004; Jackson et al. 2005), ESTER models self-consistently include the differential rotation of the radiative envelope that is driven by the baroclinic torque. They also treat self-consistently the associated meridional circulation. A brief description of these models is given in Appendix B, but we refer to the original papers of Espinosa Lara & Rieutord (2013) and Rieutord et al. (2016) for a more detailed account.

In Fig. 3 we illustrate the latitudinal variations in Γ_Ω for a 15 M_⊙ ESTER model and for a 40 M_⊙ ESTER model both taken at ZAMS, computed for the metallicity Z = 0.02 and for various values of ω that we now define as ω = Ω_eq/Ω_k (see below).

Fig. 3. Rotation-dependent Eddington parameter ΓΩ(θ) as a function of colatitude for various fractions of the Keplerian angular velocity for a 15 M⊙ (top) and a 40 M⊙ (bottom) ESTER model at ZAMS, with Z = 0.02. The grey area corresponds to supercritical rotation.

We first consider the 15 M_⊙ ESTER-model at ZAMS. Interestingly, we find that Γ_Ω(θ) first slightly decreases with increasing co-latitude before eventually vigorously increasing near equator at high angular velocity ratios. The decrease in Γ_Ω(θ) is clearly an opacity effect, which we trace back to the density decrease with θ along the stellar surface. The increase near equator at high rotation speeds is an effect of the divergence of the function f(r = 1, θ). Surprisingly, we see that Γ_Ω = 1 requires ω = 0.997, showing that the difference between the actual critical angular velocity and the Keplerian velocity is really tiny for a 15 M_⊙ ZAMS star.

To strengthen the effects of radiative acceleration, we considered the case of a 40 M_⊙ ZAMS model. Here we also see (Fig. 3 bottom) that ω must be as high as ∼0.96 for the ΩΓ limit to be reached². Technically, these latter results are not as precise as those for the 15 M_⊙ model because we approach the current limits of the ESTER code in terms of resolution, but they also point to a small difference between Ω_c and Ω_k.

From Eq. (28), at ZAMS, we find the equatorial Eddington parameter Γ_eq ≃ 0.033 at criticality for the 15 M_⊙ ESTER model and Γ_eq ≃ 0.18 for the 40 M_⊙ ESTER model. This is surprisingly low for such massive stars. To clarify this result, Fig. 4 shows the latitudinal variations in Eddington parameter Γ for both the 15 M_⊙ and 40 M_⊙ ESTER 2D models and for various angular velocity ratios. For the two models, Γ decreases with co-latitude when rotation is non-zero. The more rapid the rotation, the lower Γ_eq. The only latitudinal dependence of the Eddington parameter being on opacity, we trace back the decrease in the latter at low latitudes to the decrease in surface density (see Appendix A). We note that the opacity cannot be lower than a minimum set by pure electron scattering. While the latitudinal variations of κ are somewhat unimportant for determining the spatial location of criticality, they are crucial to the value of Ω_c/Ω_k.

Fig. 4. Eddington parameter Γ as a function of co-latitude for a 15 M⊙ (top) and a 40 M⊙ ESTER 2D model (bottom) at ZAMS and for various angular velocity ratios, with Z = 0.02.

It might be wondered, however, whether more evolved or more massive stars might have larger Γ_eq and thus a critical angular velocity that is farther from the Keplerian angular velocity. Figure 5 shows the evolution of both Γ_Ω(π/2) and Γ_eq as a function of the fractional abundance of hydrogen in the convective core X_core/X₀ for a 15 M_⊙ ESTER 2D model initially rotating at ω_i ≡ Ω_eq, i/Ω_k = 0.5, and without considering any mass loss. For the non-rotating case, evolution tends to increase Γ_eq. The increase in luminosity associated with nuclear evolution surpasses the decrease in surface opacity that is due to stellar expansion. However, when rotation is included and ω_i = 0.5, criticality is reached when X_core/X₀ ≃ 0.36. At this time, Γ_Ω diverges. While evolution proceeds and ω grows, the star flattens considerably, causing a significant drop in opacity in the equatorial region. This is clearly shown by the Γ_eq curve of Fig. 5. After a slight increase at the beginning of evolution, Γ_eq drops when criticality approaches. For this model, the equatorial opacity reduction completely dominates the effect of luminosity growth due to evolution.

Fig. 5. Γ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. The black line corresponds to the evolution of Γeq for ω = 0.

These results tend to confirm the idea put forward by Glatzel (1998), namely that the critical angular velocity is not strongly modified by the Eddington limit. We may conclude that because of the noticeable effect of rotation on opacity at the equator, the critical angular velocity is only slightly reduced compared to the Keplerian angular velocity, at least for stars with a mass lower than 40 M_⊙ at Z = 0.02.

In practice, the difference is therefore tiny enough to be neglected in view of the other uncertainties of stellar models. We therefore continue to express Ω as a fraction of the equatorial Keplerian angular velocity Ω_k to appreciate the distance to criticality, bearing in mind that this fraction is slightly smaller than the actual one.

3. Local mass-flux prescription from 1D models

We now address the second question of the paper: the dependence of mass and angular momentum losses on rotation rate. All hot stars have radiation-driven winds that become directly observable in spectral energy distributions and spectral lines as soon as they are above some luminosity threshold. For massive stars of spectral types O, B, and A, this threshold corresponds to L ∼ 10⁴ L_⊙ (Abbott 1979). Above this luminosity, massive stars show direct spectroscopic evidence of winds throughout their lifetime (UV P Cygni line profiles and optical emission lines such as Hα, see Abbott 1979 or Kudritzki & Puls 2000 and references therein). The radiative acceleration has a significant effect on the total gravity, and as shown in the previous section, reduces the critical angular velocity.

In this section we propose to derive a local mass-flux prescription that can be seen as a local equivalent of CAK original theory where the radiation-driven wind is assumed to be an isothermal stationary flow that is driven outward by photon scattering and absorption. We account for the finite cone angle of the radiating photospheric surface and for radial variations of ionisation using the results of Pauldrach et al. (1986) and Friend & Abbott (1986). Still, with or without these additional corrections, the global mass-loss rate follows a similar scaling for the wind momentum-luminosity relation (Kudritzki et al. 1995; Puls et al. 1996).

3.1. Global mass-loss rate derived from 1D CAK theory

In the 1D spherically symmetric case (i.e. without rotation), the two hydrodynamical equations needed to describe the mass-flux are the conservation of mass,

$$\begin{array}{c}\hfill \dot{M}=4\pi r^2\rho v=\mathrm{const}.,\end{array}$$(30)

and the radial momentum equation,

$$\begin{array}{c}\hfill v\frac{\partial v}{\partial r}=-\frac{1}{\rho }\frac{\partial p}{\partial r}+g+g_{\mathrm{rad}},\end{array}$$(31)

where Ṁ is the total mass-loss rate of the star, and ρ, v, p are the density, radial velocity, and gas pressure, respectively. g_rad is the radiative acceleration,

$$\begin{array}{c}\hfill g_{\mathrm{rad}}=g_{\mathrm{rad}}^{\mathrm{line}}+g_{\mathrm{e}},\end{array}$$(32)

where $g_{\mathrm{rad}}^{\mathrm{line}}$ is the line-driven acceleration, and g_e = κ_eF/c is the radiative acceleration due to Thomson scattering. κ_e is the opacity from electron scattering. Because we are interested in O, B, and A stars, both bound-free and free-free transitions are neglected (e.g. Runacres & Blomme 1994; Gayley 1995). This may not be valid for Wolf-Rayet stars, however. Finally, g is the gravitational acceleration,

$$\begin{array}{c}\hfill g=-\frac{\mathit{GM}}{r^2}·\end{array}$$(33)

Using the ideal gas equation of state, we write

$$\begin{array}{c}\hfill p=c_{\mathrm{s}}^2\rho ,\end{array}$$(34)

where c_s is the isothermal sound speed.

Because winds from hot stars are mostly line driven, the evaluation of the line-driven radiative acceleration plays a crucial role in determining the mass flux. In the Sobolev approximation (i.e. large velocity gradient approximation), considering a purely radial streaming radiation from a point-source star, the line-driven radiative acceleration can be written (CAK)

$$\begin{array}{c}\hfill g_{\text{rad}}^{\mathrm{line}}=M(t)g_{\mathrm{e}}\equiv k{\left(\frac{\partial v/\partial r}{\rho v_{\mathrm{th}}{\kappa }_{\mathrm{e}}}\right)}^{\alpha }{g}_{\mathrm{e}},\end{array}$$(35)

where M(t)=kt^−α is the CAK force multiplier. α and k are the CAK force multiplier parameters (FMPs). α can be interpreted as the ratio of the line force from optically thick lines to the total line-force, which thus decreases with decreasing effective temperature because of the increased iron group lines (e.g. Puls et al. 2000). Moreover, the quantity k is related to the fraction of the total stellar flux, which would be blocked in the photosphere if all lines were optically thick (Puls et al. 2000). t is the electron optical depth parameter and v_th is the thermal speed, usually taken as the proton thermal speed v_th ≡ (2k_BT_eff/m_H)^1/2 (e.g. Abbott 1982). However, at least in the lower part of the wind, Fe line-driving dominates. Therefore we instead take v_th ≡ (2k_BT_eff/m_Fe)^1/2 for the standard CAK formalism. We then obtain the global mass-loss rate of a non-rotating star, namely

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{CAK}}=\frac{4\pi }{{\kappa }_{\mathrm{e}}{v}_{\mathrm{th}}}{\left(\frac{k\alpha {\kappa }_{\mathrm{e}}L}{4\pi c}\right)}^{1/\alpha }{\left(\frac{1-\alpha }{\alpha }\right)}^{\frac{1-\alpha }{\alpha }}{[GM(1-{\mathrm{\Gamma }}_{\mathrm{e}})]}^{\frac{\alpha -1}{\alpha }}.\end{array}$$(36)

The global mass loss (without rotation) thus scales as

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{CAK}}\propto {[M(1-{\mathrm{\Gamma }}_{\mathrm{e}})]}^{(\alpha -1)/\alpha }{L}^{1/\alpha },\end{array}$$(37)

which is the basis of the wind momentum-luminosity relation (Kudritzki et al. 1995; Puls et al. 1996).

3.2. Finite disc and ionisation corrections

In the CAK approach, the purely radial streaming radiation leads to an electron optical depth parameter t that only depends on (dv/dr)⁻¹. This assumption neglects the finite cone angle of the radiating photospheric surface, however. Using Eq. (49) of CAK, we may rewrite t with its exact expression, leading to the modified force multiplier (Pauldrach et al. 1986), namely

$$\begin{array}{c}\hfill M(t\prime )=M(t)\frac{2}{1-{\mu }_{\ast }}{{\displaystyle \int {\mu }_{\ast }1}}_{}^{}{\left[\frac{(1-{\mu }^2)v/r+{\mu }^{2}v\prime }{v\prime }\right]}^{\alpha }\mu \mathrm{d}\mu ,\end{array}$$(38)

where v′=dv/dr, u = −R/r, and μ is the cosine of the angle between the direction of emitted radiation and the radial direction and ${\mu }_{*}=\sqrt{1-u^2}$. Evaluating the integral in Eq. (38) yields the modified force multiplier, corrected for finite cone angle, namely

$$\begin{array}{c}\hfill M(t\prime )\simeq \frac{M(t)}{u^2(1+\alpha )(1+\frac{w}{uw\prime })}[1-{(1-u^2-u\frac{w}{w\prime })}^{1+\alpha }],\end{array}$$(39)

where w = v/v_th and w′=dw/du. As in the original CAK derivation, the mass-loss rate is calculated at the critical radius r_c defined by a singularity and a regularity condition. Following Pauldrach et al. (1986), we assumed r_c to be located very close to the stellar radius, r_c ≃ R. This assumption may not be verified for rotators close to criticality, for which the fast-wind solution (with r_c ≃ R) is replaced by the so-called Ω-slow solution in the equatorial plane. This solution is characterised by an increased mass-loss rate, a slower and denser wind with a critical radius that is much farther in the wind (Curé 2004; Curé et al. 2005; Araya et al. 2017). Because the increase in mass loss associated with the Ω-slow solution was quite modest (factor ∼2 in Curé 2004), we decided to ignore it. We further assumed the velocity to follow a power law like

$$\begin{array}{c}\hfill v(u)=v_{\infty }{(1+u)}^{\beta },\end{array}$$(40)

where v_∞ is the terminal velocity of the wind and 0.7 ≲ β ≲ 1.3. The corrected force multiplier then simplifies to

$$\begin{array}{c}\hfill M(t\prime )\simeq \frac{M(t)}{1+\alpha },\end{array}$$(41)

and results in a modified prefactor for the mass-loss rate. Thus we use

$$\begin{array}{c}\hfill \dot{M}={\left(\frac{1}{1+\alpha }\right)}^{1/\alpha }{\dot{M}}_{\mathrm{CAK}}.\end{array}$$(42)

Because 0 <  α <  1, the finite-disc correction reduces the global mass-loss rate, and for typical values of α, namely between 0.4 and 0.7, the CAK mass-loss rate is multiplied by a factor ∼4/9. Additionally, the effect of radial changes in ionisation in the outward direction in the wind can be approximately taken into account by correcting the force multiplier of Eq. (41), namely multiplying it by a factor (n_e/W)^δ (Abbott 1982), where n_e is the electron density in units of 10¹¹ cm⁻³ and $W\equiv 0.5(1-\sqrt{1-u^2})$ is the radiation dilution factor. δ is then another FMP. This modification of the line-driven acceleration can be roughly accounted for by replacing α in the power exponents of Eqs. (36) and (42) with α′≡α − δ (Puls et al. 1996, 2000). Finally, we obtain the modified local mass-flux in the non-rotating case,

$$\begin{array}{cc}\hfill \dot{m}\equiv \frac{\dot{M}}{4\pi R^2}=& \left(\frac{\alpha }{v_{\mathrm{th}}c}\right){\left(\frac{k}{1+\alpha }\right)}^{1/\alpha \prime }\hfill \\ \hfill & \times {\left[\frac{c}{{\kappa }_{\mathrm{e}}(1-\alpha )}(|g|-\frac{{\kappa }_{\mathrm{e}}F}{c})\right]}^{\frac{\alpha \prime -1}{\alpha \prime }}{F}^{1/\alpha \prime },\hfill \end{array}$$(43)

where we used the radiative flux F rather than the luminosity.

Unlike the approach of MMM, we do not need to express the mass flux so that it explicitly depends on the total gravity. Rather, it now depends on gravity g, corrected for the radiative acceleration from electron scattering κ_eF/c.

3.3. Parametrisation of the FMPs

We now focus on the different FMPs α, k, and δ to estimate how they vary with the effective temperature T_eff. We assumed that δ does not significantly vary with T_eff and took δ = 0.1, a typical value for hot stars at solar metallicity (Abbott 1982). We note that δ can reach much higher values in very metal-poor stars where the wind is mostly driven by hydrogenic lines, and can even be negative under very specific conditions (Puls et al. 2000).

For α, we took fixed values at T_eff  =  10 kK, 20 kK,  30 kK, and 40 kK and imposed linear interpolation in between (J. Puls, priv. comm.), namely

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

This function is shown in Fig. 6.

Fig. 6. Adopted force multiplier parameter α as a function of the effective temperature Teff from Eq. (44) (courtesy J. Puls). The black vertical dashed lines mark the location of the imposed values for α.

Finally, we calibrated k assuming that our expression for the mass-loss rate in the non-rotating regime of Eq. (43) is equivalent to the expression of Vink et al. (2001). The Vink et al. (1999) calculations of wind models for OB stars showed that around T_eff ≃ 25 kK, the mass-loss rate Ṁ suddenly increases (towards lower T_eff) as a result of the recombination of Fe IV into Fe III, which has a stronger line acceleration in the lower part of the wind. Lamers et al. (1995) and Vink et al. (1999) suggested the existence of a second bi-stability jump, around T_eff = 10 kK, that would be caused by the recombination of Fe III into Fe II. Vink et al. (2001) did not account for this jump, however.

The Vink et al. (2001) prescription for mass loss still awaits confirmation, however. Their predictions for the size and position of the main bi-stability jump have not been confirmed by observations until today. For instance, Markova & Puls (2008) found an Ṁ jump of a factor in between 0.4 and 2.5, and more recent theoretical modelling by Petrov et al. (2016) found the bi-stability jump at $T_{\mathrm{eff}}^{\mathrm{jump}}\simeq 20\mathrm{kK}$, while Vink et al. (2001) predicted an Ṁ-jump by a factor ∼10 located at $T_{\mathrm{eff}}^{\mathrm{jump}}\simeq 25\mathrm{kK}$. Crowther et al. (2006), on the other hand, found a more gradual decrease in terminal velocity v_∞ instead (thus a more gradual increase in Ṁ). In addition, a discrepancy of a factor 2–3 also appears when the mass-loss rates of hot OB stars are compared with $T_{\mathrm{eff}}>T_{\mathrm{eff}}^{\mathrm{jump}}$ obtained with the Vink et al. (2001) models and from X-ray, UV, and IR diagnostics (e.g. Najarro et al. 2011; Sundqvist et al. 2011; Bouret et al. 2012; Cohen et al. 2013; Leutenegger et al. 2013; Hervé et al. 2013; Rauw et al. 2015). This could be due to the significant effect of small-scale inhomogeneities in the wind (e.g. Puls et al. 2008; Sundqvist & Owocki 2013; Puls et al. 2015) and/or to the outdated solar mixture used in the Vink et al. (2001) models, namely Z_⊙ ≃ 0.02 with an Anders & Grevesse (1989) mixture. The more recent solar composition with Z_⊙ ≃ 0.014 of Asplund et al. (2009) could reduce the discrepancy between predicted mass-loss rates and observations (see Sect. 4.3 for a short discussion of the effects of metallicity on mass loss). Nevertheless, the Vink et al. (2001) models are still widely used in stellar evolution codes, and we also used their recipe to calibrate k to qualitatively predict the impact of radiation-driven winds on the rotational evolution of massive stars.

We thus assumed Ṁ = Ṁ_Vink in the non-rotating case, and we calibrated k with non-rotating 1D ESTER models, that is, using mass, luminosity, and effective temperature outputs from 1D ESTER models of various masses as inputs to the Vink et al. (2001) mass loss prescription, taking v_th as the thermal velocity of Fe ions, namely v_th ≡ (2k_BT_eff/m_Fe)^1/2. Because we calibrated k using the Vink et al. (2001) mass-loss prescription, each line was considered with its appropriate v_th.

From now on, the calibrated k is written k′. We find that k′ slightly varies along the main sequence and therefore had to decide which evolution state to use for calibration. We chose to calibrate our ESTER models at ZAMS. Fitting k′ finally gives us the following semi-empirical function k′(T_eff) at Z = 0.02, defined on both sides of the effective temperature jump,

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

where Vink et al. (2001) defined

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

with ⟨ρ⟩ the characteristic wind density at 50% of the terminal velocity of the wind, given by

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

The function k′(T_eff) is shown in Fig. 7. Our values of k′ assume v_th = (2k_BT_eff/m_Fe)^1/2; other assumptions on v_th would lead to other values of k′ to remain compatible with the Vink et al. (2001) mass loss.

Fig. 7. Calibrated FMP k′ as a function of the effective temperature Teff for various 1D ZAMS models computed with ESTER for different masses with Z = 0.02. The red full line shows the corresponding fit.

4. Effects of rotation on mass and angular momentum loss

After parametrising the FMPs and expressing the local mass-flux as a function of the radiative flux as well as gravity and acceleration from free-electron scattering in the non-rotating regime, we assumed that the latter follows the same scaling laws when rotation is taken into account. We therefore ignored the changes in finite disc prefactor and set this correction to 4/9. With Eq. (43), the local mass-flux per unit surface for a rotating star reads

$$\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 (\theta ))}\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}$$(48)

This local mass-flux expression is now θ dependent and thus leads to an anisotropic stellar wind that, at first glance, would favour polar ejection due to the higher polar radiative flux³ (e.g. Owocki & Gayley 1997; Owocki et al. 1998a; Petrenz & Puls 2000; Maeder & Meynet 2000). We note that Curé (2004) took a different approach and analytically derived an equation for the mass-loss rate that accounts for rotation at the equator. This equation has a Ω-slow solution for rotators close to criticality. We now investigate the surface distribution of ṁ from the outputs of ESTER 2D models at ZAMS with Z = 0.02 and for various ω.

4.1. Latitudinal variations in mass and angular momentum loss

We computed the local mass-flux ṁ(θ) as well as the local angular momentum flux

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

with ESTER 2D models and prescription (48).

Figure 8 shows the local mass-flux ṁ and local angular momentum flux $\dot{\ell }$ as a function of colatitude for ESTER 2D models of a 15 M_⊙ star with Z = 0.02, at ZAMS and for various angular velocity ratios ω = Ω_eq/Ω_k. When the star is far from the Ω-limit, mass loss is favoured in polar regions and thus decreases towards the equator. However, for a sufficiently high angular velocity ratio (ω ≳ 0.85 in Fig. 8), there is a colatitude θ_jump where $T_{\text{eff}}({\theta }_{\mathrm{jump}})=T_{\text{eff}}^{\text{jump}}$ and the local bi-stability limit is crossed. In that case, the mass flux is enhanced between θ_jump and the equator and the star is in a two-wind regime (TWR), otherwise it is in a single-wind regime (SWR). This local bi-stability jump therefore strongly modifies the distribution of the mass flux with colatitude: while mass loss is favoured in polar regions in the SWR, it is far stronger in equatorial regions in the TWR. Similarly, while the angular momentum flux is maximum at some intermediate colatitude in the SWR, it is strongly favoured in the equatorial regions in the TWR. The idea of an enhanced equatorial mass-flux that is due to both gravity darkening and the local bi-stability limit has also been discussed in the past, for instance, by Zickgraf et al. (1986, 1989), Lamers & Pauldrach (1991), Owocki & Gayley (1997), Owocki et al. (1998b), and Pelupessy et al. (2000).

Fig. 8. Variation in surface mass flux ṁ (top) and surface angular momentum flux $\dot{\ell }$ (bottom) as a function of colatitude θ for a 15 M⊙ ESTER 2D-model at ZAMS with Z = 0.02 and various angular velocity ratios ω. The main bi-stability limit is reached near the equator for ω ≳ 0.85, and induces a strong mass-flux and angular momentum flux for colatitudes in between θjump and the equator.

This change in the latitudinal distribution of the angular momentum flux is particularly important for stellar evolution. In the SWR, polar-dominated mass loss allows rapidly rotating massive stars to lose mass during the MS without losing much angular momentum, hence keeping a rapid rotation throughout their evolution. In the TWR, however, mass loss becomes highly dominated by the equatorial regions and the star loses far more angular momentum. This enhanced loss of angular momentum in the TWR could therefore induce a drop in ω during stellar evolution. This phenomenon will be discussed in the follow-up paper and is not to be confused with the bi-stability braking introduced by Vink et al. (2010), which is purely one-dimensional and corresponds to the global transition between the hot and cold side of the bi-stability jump. We note that a star need not be close to Keplerian rotation to reach the local bi-stability limit. A rotating star that has an equatorial effective temperature that is only slightly higher than the temperature of the jump can reach the TWR with a small increase of ω.

4.2. Effects of rotation on the global mass and angular momentum loss rates

We now compute the global mass and angular momentum loss rates by integrating ṁ(θ) and $\dot{\ell }(\theta )$ over the distorted stellar surface as follows:

$$\begin{array}{c}\hfill \dot{M}=2\pi {\displaystyle \int \dot{m}(\theta )R^2(\theta )\sqrt{1+\frac{R_{\theta }^2}{R^2(\theta )}}sin\theta \mathrm{d}\theta ,}\end{array}$$(50)

$$\begin{array}{c}\hfill \dot{\mathcal{L}}=2\pi {\displaystyle \int \dot{\ell }(\theta )R^2(\theta )\sqrt{1+\frac{R_{\theta }^2}{R^2(\theta )}}sin\theta \mathrm{d}\theta ,}\end{array}$$(51)

where R(θ) is the θ-dependent radius of the star. The area element at the stellar surface is

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

where R_θ = ∂R/∂θ (Rieutord et al. 2016).

The global mass-loss rate Ṁ, the critical angular velocity ratio ω_c = Ω_eq/Ω_c as given by the ω-model (Eq. (28)), the ratio of equatorial angular velocity to Keplerian angular velocity ω = Ω_eq/Ω_k, and the angular momentum loss timescale ${\mathrm{T}}_L=\mathcal{L}/\dot{\mathcal{L}}$ are reported in Table 1 for a 15 M_⊙ star ESTER 2D model at ZAMS and at Z = 0.02. For this stellar model in the SWR, we find the global mass-loss rate to slightly decrease for increasing ω, for instance, Ṁ(ω = 0.6) / Ṁ(ω = 0) ≃ 0.87 (see Fig. 9 top). Similar results have been obtained by Müller & Vink (2014).

Fig. 9. Variation in mass loss rate Ṁ (in M⊙ yr−1, top) and the angular momentum loss timescale Tℒ (in yr, bottom) as a function of the angular velocity ratio ω for a 15 M⊙ star at ZAMS with Z = 0.02. The dashed lines show the same as the solid line, but with an FMP α′ that has been reduced by 1% to show the sensitivity of Ṁ and Tℒ to FMP variations.

With increasing ω, the total angular momentum of the star ℒ increases, and even though the global mass-loss rate Ṁ decreases in the SWR, the global loss of angular momentum $\dot{\mathcal{L}}$ also increases in this regime. This is simply because $\dot{\mathcal{L}}$ increases for increasing ω. It is even more interesting, however, that the timescale of angular momentum loss $T_{\mathcal{L}}=\mathcal{L}/\dot{\mathcal{L}}$ is approximately independent of the degree of criticality ω in the SWR (see Fig. 9, bottom).

On the other hand, the TWR is characterised by a strong increase in global mass and angular momentum loss rates. In this regime, T_ℒ rapidly decreases as ω approaches unity. Both the strong increase in Ṁ and decrease in T_ℒ can be explained by the increasing stellar surface fraction where the effective temperature is lower than $T_{\mathrm{eff}}^{\mathrm{jump}}$ as ω increases (see Fig. 8).

That Ṁ only gradually increases with increasing ω in the TWR is a result specific to 2D models. In 1D models, the bi-stability jump is accounted for with a stronger global mass-loss rate if the mean effective temperature of the star is lower than $T_{\mathrm{eff}}^{\mathrm{jump}}\sim 22.5$–25 kK, according to Vink et al. (2001). In the present work however, 2D models reach the bi-stability limit if the local effective temperature somewhere on the stellar surface is lower than $T_{\mathrm{eff}}^{\mathrm{jump}}$. This difference has two consequences. Firstly, 2D models can reach the bi-stability limit even with an average effective temperatures higher than $T_{\mathrm{eff}}^{\mathrm{jump}}$. In Fig. 10, we illustrate the global mass-loss rate for a variety of angular velocity ratios at ZAMS for a 15 M_⊙ and a 10 M_⊙ ESTER model at ZAMS and with Z = 0.02, against the corresponding surface-averaged effective temperature ${\overline{T}}_{\mathrm{eff}}$ of the model. In these models, ${\overline{T}}_{\mathrm{eff}}$ is greater than $T_{\mathrm{eff}}^{\mathrm{jump}}\simeq 22.8\mathrm{kK}$ for all ω. Thus, equivalent 1D models would just have ignored the bi-stability jump.

Fig. 10. Global mass-loss rate Ṁ (in M⊙ yr−1) for a 15 M⊙ (red) and a 10 M⊙ (green) ESTER model against the corresponding mean effective temperature ${\overline{T}}_{\mathrm{eff}}$, at ZAMS with Z = 0.02 and for ω ∈ [0, 1[. Arrows indicate the direction of increasing ω.

Secondly, in 2D models the surface fraction where the effective temperature is lower than $T_{\mathrm{eff}}^{\mathrm{jump}}$ monotonically increases with increasing ω. This results in a gradual variation in Ṁ (and T_ℒ) with ${\overline{T}}_{\mathrm{eff}}$ in the TWR (see Fig. 10). Hence, in rotating stars the bi-stability jump does not induce a discontinuity of the global mass-loss rate (but it induces a discontinuity of the local mass-flux, see Fig. 8) as the ω parameter increases (and therefore as the mean effective temperature decreases). The discontinuity occurs only on the derivative of the function Ṁ(T_eff). This is further discussed in the follow-up paper.

These points show that even though the bi-stability jump might eventually be confirmed observationally (although its location in terms of mean effective temperature would be scattered, see Fig. 10), a full 2D spectral analysis is required to verify both qualitative and quantitative features of radiation-driven winds from rapidly rotating massive stars (e.g. Petrenz & Puls 1996). As a first step, a way around this full analysis would be to select stars with a small vsini to select either slowly rotating stars or stars that are viewed pole-on. Obviously, the precise determination of vsini for hot massive stars is a challenge in particular because rotational effects are mixed with other line-broadening effects such as the so-called macro-turbulence (e.g. Simón-Díaz & Herrero 2007).

4.3. Metallicity effect

Before we conclude this paper, a few words on low-metallicity stars are in order. Metallicity is indeed known to have an important effect on radiatively driven winds because metallic lines, which significantly contribute to opacity, weaken and eventually disappear. As a consequence, the FMPs, such as α or k, are expected to decrease with a decreasing Z (Kudritzki et al. 1987; Puls et al. 2000, 2008). Moreover, mass loss is very sensitive to the value of the FMPs. In particular, a small decrease in α′ leads to a significant decrease in Ṁ (see Fig. 9). In addition, a low metallicity causes stars to be more compact and therefore have a higher effective temperature (Maeder 2009). This effect may compensate (partly?) for the loss of opacity on the wind acceleration. All in all, because the FMPs have a significant influence on mass-loss calculations and because they are ill-known at metallicities much lower than solar, we do not venture any prediction on the behaviour of mass flux at low Z. We leave this question to future investigations.

5. Discussion and conclusions

We investigated two questions that are a prerequisite to the study of the evolution of massive rapidly rotating stars: (i) What is the critical angular velocity of a star when radiative acceleration is significant in its atmosphere? (ii) How do the mass and angular momentum loss rates depend on the stellar rotation rate?

To the first question, we answer that the critical angular velocity is very close to the classical Keplerian angular velocity at the equator, at least for stars with masses lower than 40 M_⊙ (and for Z = 0.02). The role of radiative acceleration turns out to be rather limited because of the combination of a reduced opacity and reduced flux at the equator. The reduction of the flux, the so-called gravity darkening, is less than was predicted by the von Zeipel model. This latter point is the main difference between this study and the pioneering investigations of Maeder (1999) and Maeder & Meynet (2000). ESTER 2D models indeed show that the flux is almost anti-parallel to gravity in the stellar radiative envelope (Espinosa Lara & Rieutord 2011). To a very good approximation, we can therefore write F = f(r, θ)g_eff, which is the base of the ω-model (Espinosa Lara & Rieutord 2011; Rieutord et al. 2016). We showed that the ω-model remains close to full 2D ESTER up to rotation as high as 90% of the critical rotation. When equatorial rotation approaches Keplerian rotation, f(r, θ) diverges at the equator, while in the von Zeipel model it remains finite. This means that the effective temperature decreases more slowly at the stellar equator than what is predicted with the von Zeipel recipe. f(r, θ) is also a monotonically increasing function of co-latitude. Its maximum is therefore reached at the equator, hence it turns out that the total acceleration g_tot = g_rad + g_eff vanishes first at the equator, when rotation is increased. Unlike the von Zeipel approximation, ESTER 2D models never predict that the radiative flux vanishes at the equator. Critical rotation, defined as the rotation required for g_tot to vanish somewhere at the surface, is therefore always reached before the equatorial rotation reaches Keplerian rotation. This point has been made by Maeder & Meynet (2000). However, 2D models hold that this difference is tiny. Considering a massive star of 40 M_⊙, we therefore find that criticality, g_tot = 0 at the equator, is reached at ω ∼ 0.96 and even at 0.997 for a 15 M_⊙ star. This tiny difference can be understood because gravity darkening in the ω-model is weaker than in the von Zeipel model and because the effect of rotation on opacity leads to a strong decrease in standard Eddington parameter towards the equator. To return to the debate between Langer (1997, 1998) and Glatzel (1998), our results support the latter concerning the influence of the Eddington limit on the value of critical rotation: this influence is quite small and never exceeds 4% as far as we could test (i.e. M ≤ 40 M_⊙, Z = 0.02). The fact that only a small equatorial region becomes unbound at criticality may lead to mechanical mass loss. This will be discussed in a forthcoming work.

To address the second question, we first devised a prescription for the surface density of the mass flux based on current knowledge of radiatively driven winds. The derivation of this local mass flux was based on the approaches of Castor et al. (1975) and Pauldrach et al. (1990), but force multiplier parameters were adjusted to match the widely used prescriptions of Vink et al. (2001) for Ṁ in the range T_eff ∈ [10, 50] kK. This prescription led to a discontinuity in the mass flux when T_eff drops below 22.5–25 kK. Because the surface effective temperature of rotating stars can span a wide range of values from poles to equator, it easily happens that the discontinuity occurs at some latitude of the star. In this case, the stellar wind experiences two regimes, one centred on the poles, the other around the equator. We have shown that if the star experiences a single-wind regime (no latitude of discontinuity), the maximum extraction of angular momentum occurs at mid-latitude, while the mass flux is maximum at the poles. However, if the jump in mass flux occurs at some latitude, then both mass loss and angular momentum loss are maximum in equatorial regions. Interestingly, these two regimes are expected to affect not only the evolution of the stellar rotation rate, but also the internal rotational mixing because the applied torque is different in both intensity and location.

Before we conclude, we wish to caution about one important simplification of ESTER 2D models. The current ESTER models indeed assume that no mass flux leaves the photosphere and a zero normal velocity is imposed at the surface of the star. Moreover, we assume that as in 1D models, the surface layers are vertically in hydrostatic equilibrium. All these approximations are acceptable for determining the bulk structure of the star, but are likely too rough to describe the surface layers of a wind-emitting massive star. In particular, the values of the surface opacity, which is important for determining the radiative acceleration, may be modified when a better coupling between the wind and the star is introduced. With such a new 2D model of the wind launch region, the concept and conditions of critical angular velocity will have to be revisited. With the current models, predictions are therefore indicative: they are reliable for intermediate-mass stars (lower than 10 M_⊙), but their realism and their reliability decrease with increasing mass. Beyond 40 M_⊙, new models are probably mandatory to obtain a sensible description of the mass-loss phenomenon with rotation.

Finally, on the observational side, we remark that rotation makes verifying the existence of the jump in the relation Ṁ(T_eff) more difficult. This verification would be possible if we could select stars whose rotation axis is aligned with the line of sight. In that case, we would be sure to face the same (polar) wind regime. If no selection can be made, the random orientation of the rotation axis means that the observed winds are sourced by an unconstrained range of T_eff, implying that any discontinuity in the Ṁ(T_eff) relation is smoothed out, unless we can reproduce the observed star with a complete 2D wind+star model. In a follow-up paper (Gagnier et al. 2019), we apply these results to study the evolution of rotation in early-type fast-rotating stars and address the question, among others, how a wind can prevent a massive star from reaching the critical rotation.

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: Opacity dependence on effective temperature

In this appendix, we show that gravity darkening at the surface of rotating stars may lead to a decrease in opacity κ towards the equator. To do this, we assumed that the opacity at the stellar surface follows Kramer’s opacity law, namely,

$$\begin{array}{c}\hfill \kappa (\theta )\propto {\rho }_{\mathrm{s}}(\theta )T_{\mathrm{eff}}{(\theta )}^{-7/2},\end{array}$$(A.1)

where ρ_s(θ) is the local surface density. This approximation seems to be rather well verified at the surface of rapidly rotating ESTER 2D models. We recall that the pressure at the surface is

$$\begin{array}{c}\hfill P_{\mathrm{s}}(\theta )={\tau }_{\mathrm{s}}\frac{g_{\mathrm{eff}}(\theta )}{\kappa (\theta )}\propto \frac{g_{\mathrm{eff}}(\theta )}{\kappa (\theta )},\end{array}$$(A.2)

where τ_s ≃ 2/3 is the Rosseland mean optical depth at the photosphere. Assuming the power law $T_{\text{eff}}\propto g_{\text{eff}}^{\beta }$ and the ideal gas equation of state, the previous expression leads to

$$\begin{array}{c}\hfill \kappa (\theta )\propto T_{\mathrm{eff}}{(\theta )}^{-9/4+1/(2\beta )}\mathrm{with}{\rho }_{\mathrm{s}}(\theta )\propto T_{\mathrm{eff}}{(\theta )}^{5/4+1/(2\beta )}.\end{array}$$(A.3)

The value β = 0.25 given by von Zeipel’s law implies that the opacity increases towards the equator, that is, with decreasing effective temperature. When β <  2/9 ≃ 0.222, however, this simple model shows that the surface opacity decreases towards the equator. Because β decreases with rotation, β <  2/9 corresponds to a surface flattening ϵ ≳ 0.08 (or to an angular velocity ratio ω ≳ 0.4, according to ESTER 2D models). This scaling relation is only approximate. Still, it shows that a weaker gravity darkening than that of von Zeipel may have a strong effect on the latitudinal variations in surface density, thus on opacity at the surface of rotating stars. In other words, in some cases, rotation may induce a decrease in opacity towards the equator because of a corresponding reduced density in these regions.

Appendix B: Short presentation of ESTER models

The ESTER code computes the steady state of an isolated rotating star, including the large-scale flows driven by the baroclinicity of the radiative regions. It solves in two dimensions (assuming axisymmetry) the steady equations of stellar structure with fluid flows, namely the Poisson equation,

$$\begin{array}{c}\hfill \mathrm{\Delta }\varphi =4\pi G\rho ,\end{array}$$(B.1)

where ϕ is is the gravitational potential; the continuity equation,

$$\begin{array}{c}\hfill \mathbf{\nabla }·\rho \mathit{v}=0,\end{array}$$(B.2)

the momentum equation,

$$\begin{array}{c}\hfill \rho \mathit{v}·\mathbf{\nabla }\mathit{v}=-\mathbf{\nabla }P-\rho \mathbf{\nabla }\varphi +{\mathit{F}}_{\mathrm{visc}},\end{array}$$(B.3)

where F_visc is the viscous force; and the heat balance equation,

$$\begin{array}{c}\hfill \rho T\mathit{v}·\mathbf{\nabla }s=\mathbf{\nabla }·(\chi \mathbf{\nabla }T)+{\epsilon }_{\ast }\mathrm{in}\mathrm{radiative}\mathrm{envelopes}\end{array}$$(B.4)

and

$$\begin{array}{c}\hfill {\partial }_{r}s=0\mathrm{in}\mathrm{convective}\mathrm{cores}.\end{array}$$(B.5)

This last equation assumes an efficient convection in convective cores, as can be shown with the mixing-length model.

These equations are completed by boundary conditions that require that (i) the gravitational potential ϕ vanishes at infinity, (ii) velocity fields meet stress-free conditions at the stellar surface, (iii) that the surface radiates like a local black body, and (iv) the surface is defined by the place where the pressure P equals g_eff/κ. Usual notations have been used: G is the gravitational constant, v the velocity field, s the entropy, and ε_* the energy produced by nuclear reactions per unit mass.

Regarding the micro-physics, opacity and the equation of state are given by the OPAL tables using the GN93 mixture (Grevesse & Noels 1993). It might be argued that the use of the GN93 mixture is questionable considering that a newer solar chemical composition have been determined (e.g. Asplund et al. 2009; Przybilla et al. 2009), but it is sufficient because this newer composition is not so different from the solar mixture used in Vink et al. (2001) (namely Anders & Grevesse 1989). The diffusive transport of momentum is ensured by a vanishingly low viscosity, implying that no heat is advected by meridional circulation (this is the zero Prandtl number limit). However, differential rotation resulting from the baroclinic torque is computed as well as the associated meridional circulation. Nuclear energy generation is described by an analytical formula including the pp- and CNO cycles. A more detailed description can be found in Espinosa Lara & Rieutord (2013) and Rieutord et al. (2016).

All Tables

All Figures

	Fig. 1. Scaled total gravity, effective gravity, and radiative acceleration at the equator as a function of ω with Γeq = 0.5. The total gravity at the equator has two zeros; the first root corresponds to ΓΩ(π/2)=1, and the second to geff(π/2)=0.
In the text

Fig. 2. Relative difference between the radiative flux of the ω-model and the ESTER model as a function of co-latitude for a 15 M⊙ (top) and a 40 M⊙ (bottom) ZAMS-star with various angular velocity ratios. The black line corresponds to an evolved 15 M⊙ ESTER 2D model with ω = 0.9 and a fractional abundance of hydrogen in the convective core Xcore/X0 = 0.5. X0 is the initial hydrogen mass fraction at ZAMS. The minimum of each curve corresponds to a sign change of FESTER − Fω.

In the text

	Fig. 3. Rotation-dependent Eddington parameter ΓΩ(θ) as a function of colatitude for various fractions of the Keplerian angular velocity for a 15 M⊙ (top) and a 40 M⊙ (bottom) ESTER model at ZAMS, with Z = 0.02. The grey area corresponds to supercritical rotation.
In the text

	Fig. 4. Eddington parameter Γ as a function of co-latitude for a 15 M⊙ (top) and a 40 M⊙ ESTER 2D model (bottom) at ZAMS and for various angular velocity ratios, with Z = 0.02.
In the text

	Fig. 5. Γ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. The black line corresponds to the evolution of Γeq for ω = 0.
In the text

	Fig. 6. Adopted force multiplier parameter α as a function of the effective temperature Teff from Eq. (44) (courtesy J. Puls). The black vertical dashed lines mark the location of the imposed values for α.
In the text

	Fig. 7. Calibrated FMP k′ as a function of the effective temperature Teff for various 1D ZAMS models computed with ESTER for different masses with Z = 0.02. The red full line shows the corresponding fit.
In the text

	Fig. 8. Variation in surface mass flux ṁ (top) and surface angular momentum flux $\dot{\ell }$ (bottom) as a function of colatitude θ for a 15 M⊙ ESTER 2D-model at ZAMS with Z = 0.02 and various angular velocity ratios ω. The main bi-stability limit is reached near the equator for ω ≳ 0.85, and induces a strong mass-flux and angular momentum flux for colatitudes in between θjump and the equator.
In the text

	Fig. 9. Variation in mass loss rate Ṁ (in M⊙ yr−1, top) and the angular momentum loss timescale Tℒ (in yr, bottom) as a function of the angular velocity ratio ω for a 15 M⊙ star at ZAMS with Z = 0.02. The dashed lines show the same as the solid line, but with an FMP α′ that has been reduced by 1% to show the sensitivity of Ṁ and Tℒ to FMP variations.
In the text

	Fig. 10. Global mass-loss rate Ṁ (in M⊙ yr−1) for a 15 M⊙ (red) and a 10 M⊙ (green) ESTER model against the corresponding mean effective temperature ${\overline{T}}_{\mathrm{eff}}$, at ZAMS with Z = 0.02 and for ω ∈ [0, 1[. Arrows indicate the direction of increasing ω.
In the text