Issue 
A&A
Volume 625, May 2019



Article Number  A88  
Number of page(s)  14  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201834599  
Published online  17 May 2019 
Critical angular velocity and anisotropic mass loss of rotating stars with radiationdriven winds
^{1}
IRAP, Université de Toulouse, CNRS, UPS, CNES, 14, avenue Édouard Belin, 31400 Toulouse, France
email: damien.gagnier@irap.omp.eu
^{2}
Department of Astronomy, University of Geneva, Chemin des Maillettes 51, 1290 Versoix, Switzerland
^{3}
Space Research Group, University of Alcalá, 28871 Alcalá de Henares, Spain
Received:
7
November
2018
Accepted:
1
April
2019
Context. The understanding of the evolution of earlytype stars is tightly related to that of the effects of rapid rotation. For massive stars, rapid rotation combines with their strong radiationdriven wind.
Aims. The aim of this paper is to investigate 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 depend on the rotation rate?
Methods. To investigate fast rotation, which makes stars oblate, we used the 2D ESTER models and a simplified approach, the ωmodel, which gives the latitudinal dependence of the radiative flux in a centrifugally flattened radiative envelope.
Results. We find that radiative acceleration only mildly influences the critical angular velocity, at least for stars with masses lower than 40 M_{⊙}. For instance, a 15 M_{⊙} star on the zeroage main sequence would reach criticality at a rotation rate equal to 0.997 the Keplerian equatorial rotation rate. We explain this mild reduction of the critical angular velocity compared to the classical Keplerian angular velocity by the combined effects of gravity darkening and a reduced equatorial opacity that is due to the centrifugal acceleration. To answer the second question, we first devised a model of the local surface mass flux, which we calibrated with previously developed 1D models. The discontinuity (the socalled bistability jump) included in the Ṁ − T_{eff} relation of 1D models means that the mass flux of a fastrotating star is controlled by either a single wind or a twowind regime. Mass and angular momentum losses are strong around the equator if the star is in the twowind regime. We also show that the difficulty of selecting massive stars that are viewed poleon makes detecting the discontinuity in the relation between mass loss and effective temperature also quite challenging.
Key words: stars: rotation / stars: massloss / stars: earlytype
© D. Gagnier et al. 2019
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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 largescale 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 smallscale turbulence and largescale circulation induced by rotation in radiative zones is inserted in 1D evolution models either as an advectiondiffusion 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 (socalled rotationinduced 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). Earlytype stars are, however, often considered to be fast rotators. The hypotheses and approximations of current prescriptions are therefore no longer valid for such stars. Betype stars, for instance, are well known to be fast rotators close to the breakup 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 breakup 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, earlytype stars may also be very luminous and therefore have high radiation pressure at their surface. The induced radiationdriven 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 selfconsistent 2D models of rapidly rotating earlytype 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 rotationinduced 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 fastrotating stars. To this end, we first focus on deriving a local massflux prescription based on the 1D CAK (Castor et al. 1975) and mCAK (Pauldrach et al. 1986) theories for nonrotating 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,
where effective gravity (or gravito–centrifugal acceleration) g_{eff} is supplemented by radiative acceleration,
where κ is the fluxweighted opacity^{1}, which we approximate with the total Rosseland mean opacity, F the radiative flux, and c the speed of light.
The socalled ΩΓ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
which is the breakup 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
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) reinvestigated 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 rotationdependent 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
where χ = 4acT^{3}/(3κρ) is the radiative conductivity. Additionally, assuming solidbody rotation, Zahn (1992) obtained
where
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
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
We note that Eq. (8) assumes solidbody rotation, while Eq. (9) corresponds to the case of slow rotation (Zahn 1992). Maeder (1999) noted that ζ(θ)∼10^{−2} 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
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 earlytype star can be written as
where f(r, θ) is some function of the position to be determined. In this assumption, F and g_{eff} are assumed to be antiparallel. 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 firstorder equation of the flux can be completed by the boundary condition
The equation for f(r, θ) can then be solved analytically (see Espinosa Lara & Rieutord 2011; Rieutord 2016, for details), with the following result:
where ψ(r, θ) is obtained by solving
In this equation r has been scaled with the equatorial radius R_{eq} and
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,
so that the equatorial radiative flux reads
In the slow rotation limit, Eq. (17) can be written
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
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
From this expression, we define the rotationdependent Eddington parameter Γ_{Ω}(θ) as the ratio of the actual flux F(θ) obtained with the ωmodel, and the limiting flux, namely
Using Eq. (21) we can rewrite Eq. (1) as
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 colatitude, 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
where
and
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 as
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})^{−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 g_{eff}(π/2)=0. 
For subcritical 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 .
In line with the ωmodel and the maximum of Γ_{Ω} at equator, the condition giving the critical angular velocity Ω_{c} is
or equivalently,
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 . 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 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 fastrotating stars all show that β < 1/4 when the surface flux distribution is assumed to vary as . 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 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 antiparallel 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 zeroage main sequence (ZAMS) models of 15 M_{⊙} and 40 M_{⊙}, we computed the relative difference between the fluxes of the ESTER and ωmodel, namely
as a function of colatitude. 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 colatitude for a 15 M_{⊙} (top) and a 40 M_{⊙} (bottom) ZAMSstar 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 X_{core}/X_{0} = 0.5. X_{0} is the initial hydrogen mass fraction at ZAMS. The minimum of each curve corresponds to a sign change of F_{ESTER} − 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 followup 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 midMS 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 earlytype star. Compared to previous attempts of making stellar models in two dimensions (e.g. Roxburgh 2004; Jackson et al. 2005), ESTER models selfconsistently include the differential rotation of the radiative envelope that is driven by the baroclinic torque. They also treat selfconsistently 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. Rotationdependent 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_{⊙} ESTERmodel at ZAMS. Interestingly, we find that Γ_{Ω}(θ) first slightly decreases with increasing colatitude 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^{2}. 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 colatitude when rotation is nonzero. 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 colatitude 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_{0} for a 15 M_{⊙} ESTER 2D model initially rotating at ω_{i} ≡ Ω_{eq, i}/Ω_{k} = 0.5, and without considering any mass loss. For the nonrotating 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} ≃ 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 massflux 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 radiationdriven 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^{4} 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 massflux prescription that can be seen as a local equivalent of CAK original theory where the radiationdriven 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 massloss rate follows a similar scaling for the wind momentumluminosity relation (Kudritzki et al. 1995; Puls et al. 1996).
3.1. Global massloss rate derived from 1D CAK theory
In the 1D spherically symmetric case (i.e. without rotation), the two hydrodynamical equations needed to describe the massflux are the conservation of mass,
and the radial momentum equation,
where Ṁ is the total massloss rate of the star, and ρ, v, p are the density, radial velocity, and gas pressure, respectively. g_{rad} is the radiative acceleration,
where is the linedriven acceleration, and g_{e} = κ_{e}F/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 boundfree and freefree transitions are neglected (e.g. Runacres & Blomme 1994; Gayley 1995). This may not be valid for WolfRayet stars, however. Finally, g is the gravitational acceleration,
Using the ideal gas equation of state, we write
where c_{s} is the isothermal sound speed.
Because winds from hot stars are mostly line driven, the evaluation of the linedriven 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 pointsource star, the linedriven radiative acceleration can be written (CAK)
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 lineforce, 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_{B}T_{eff}/m_{H})^{1/2} (e.g. Abbott 1982). However, at least in the lower part of the wind, Fe linedriving dominates. Therefore we instead take v_{th} ≡ (2k_{B}T_{eff}/m_{Fe})^{1/2} for the standard CAK formalism. We then obtain the global massloss rate of a nonrotating star, namely
The global mass loss (without rotation) thus scales as
which is the basis of the wind momentumluminosity 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)^{−1}. 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
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 . Evaluating the integral in Eq. (38) yields the modified force multiplier, corrected for finite cone angle, namely
where w = v/v_{th} and w′=dw/du. As in the original CAK derivation, the massloss 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 fastwind solution (with r_{c} ≃ R) is replaced by the socalled Ωslow solution in the equatorial plane. This solution is characterised by an increased massloss 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
where v_{∞} is the terminal velocity of the wind and 0.7 ≲ β ≲ 1.3. The corrected force multiplier then simplifies to
and results in a modified prefactor for the massloss rate. Thus we use
Because 0 < α < 1, the finitedisc correction reduces the global massloss rate, and for typical values of α, namely between 0.4 and 0.7, the CAK massloss 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^{11} cm^{−3} and is the radiation dilution factor. δ is then another FMP. This modification of the linedriven 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 massflux in the nonrotating case,
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 κ_{e}F/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 metalpoor 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
This function is shown in Fig. 6.
Fig. 6. Adopted force multiplier parameter α as a function of the effective temperature T_{eff} 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 massloss rate in the nonrotating 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 massloss 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 bistability 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 bistability 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 bistability jump at , while Vink et al. (2001) predicted an Ṁjump by a factor ∼10 located at . 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 massloss rates of hot OB stars are compared with obtained with the Vink et al. (2001) models and from Xray, 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 smallscale 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 massloss 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 radiationdriven winds on the rotational evolution of massive stars.
We thus assumed Ṁ = Ṁ_{Vink} in the nonrotating case, and we calibrated k with nonrotating 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_{B}T_{eff}/m_{Fe})^{1/2}. Because we calibrated k using the Vink et al. (2001) massloss 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 semiempirical function k′(T_{eff}) at Z = 0.02, defined on both sides of the effective temperature jump,
where Vink et al. (2001) defined
with ⟨ρ⟩ the characteristic wind density at 50% of the terminal velocity of the wind, given by
The function k′(T_{eff}) is shown in Fig. 7. Our values of k′ assume v_{th} = (2k_{B}T_{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 T_{eff} 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 massflux as a function of the radiative flux as well as gravity and acceleration from freeelectron scattering in the nonrotating 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 massflux per unit surface for a rotating star reads
This local massflux 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^{3} (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 massloss 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 massflux ṁ(θ) as well as the local angular momentum flux
with ESTER 2D models and prescription (48).
Figure 8 shows the local massflux ṁ and local angular momentum flux 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 and the local bistability limit is crossed. In that case, the mass flux is enhanced between θ_{jump} and the equator and the star is in a twowind regime (TWR), otherwise it is in a singlewind regime (SWR). This local bistability 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 massflux that is due to both gravity darkening and the local bistability 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 (bottom) as a function of colatitude θ for a 15 M_{⊙} ESTER 2Dmodel at ZAMS with Z = 0.02 and various angular velocity ratios ω. The main bistability limit is reached near the equator for ω ≳ 0.85, and induces a strong massflux 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, polardominated 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 followup paper and is not to be confused with the bistability braking introduced by Vink et al. (2010), which is purely onedimensional and corresponds to the global transition between the hot and cold side of the bistability jump. We note that a star need not be close to Keplerian rotation to reach the local bistability 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 over the distorted stellar surface as follows:
where R(θ) is the θdependent radius of the star. The area element at the stellar surface is
where R_{θ} = ∂R/∂θ (Rieutord et al. 2016).
The global massloss 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 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 massloss 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).
Summary of the main results for ESTER 2D models of a 15 M_{⊙} star with Z = 0.02 at ZAMS.
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 massloss rate Ṁ decreases in the SWR, the global loss of angular momentum also increases in this regime. This is simply because increases for increasing ω. It is even more interesting, however, that the timescale of angular momentum loss 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 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 bistability jump is accounted for with a stronger global massloss rate if the mean effective temperature of the star is lower than –25 kK, according to Vink et al. (2001). In the present work however, 2D models reach the bistability limit if the local effective temperature somewhere on the stellar surface is lower than . This difference has two consequences. Firstly, 2D models can reach the bistability limit even with an average effective temperatures higher than . In Fig. 10, we illustrate the global massloss 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 surfaceaveraged effective temperature of the model. In these models, is greater than for all ω. Thus, equivalent 1D models would just have ignored the bistability jump.
Fig. 10. Global massloss rate Ṁ (in M_{⊙} yr^{−1}) for a 15 M_{⊙} (red) and a 10 M_{⊙} (green) ESTER model against the corresponding mean effective temperature , 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 monotonically increases with increasing ω. This results in a gradual variation in Ṁ (and T_{ℒ}) with in the TWR (see Fig. 10). Hence, in rotating stars the bistability jump does not induce a discontinuity of the global massloss rate (but it induces a discontinuity of the local massflux, 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 followup paper.
These points show that even though the bistability 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 radiationdriven 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 poleon. Obviously, the precise determination of vsini for hot massive stars is a challenge in particular because rotational effects are mixed with other linebroadening effects such as the socalled macroturbulence (e.g. SimónDíaz & Herrero 2007).
4.3. Metallicity effect
Before we conclude this paper, a few words on lowmetallicity 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 massloss calculations and because they are illknown 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 socalled 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 antiparallel 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 colatitude. 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 singlewind regime (no latitude of discontinuity), the maximum extraction of angular momentum occurs at midlatitude, 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 windemitting 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 intermediatemass 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 massloss 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 followup paper (Gagnier et al. 2019), we apply these results to study the evolution of rotation in earlytype fastrotating stars and address the question, among others, how a wind can prevent a massive star from reaching the critical rotation.
ω has to be slightly lower than 0.96 so that Γ_{Ω} is exactly unity at equator. At ω = 0.96, the star is already supercritical at the equator (grey area in Fig. 3).
Our approach implicitly assumes the presence of a weak, polewards directed component of the radiation force. Such a nonradial component is the result of the decreasing radial velocity towards the equator, and is essential for inhibiting a flow that otherwise would be directed towards the equator. Within our approach, however, this component can be neglected when estimating the thetadependence of ṁ. For details, see Owocki et al. (1998a), for example.
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 2017P0107). M. Rieutord acknowledges the strong support of the French Agence Nationale de la Recherche (ANR), under grant ESRR (ANR16CE31000701), 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 ESP201788436R.
References
 Abbott, D. C. 1979, in Mass Loss and Evolution of OType Stars, eds. P. S. Conti, & C. W. H. De Loore, IAU Symp., 83, 237 [NASA ADS] [CrossRef] [Google Scholar]
 Abbott, D. C. 1982, ApJ, 259, 282 [NASA ADS] [CrossRef] [Google Scholar]
 Amard, L., Palacios, A., Charbonnel, C., Gallet, F., & Bouvier, J. 2016, A&A, 587, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197 [Google Scholar]
 Araya, I., Jones, C. E., Curé, M., et al. 2017, ApJ, 846, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Bastian, N., CabreraZiri, I., Niederhofer, F., et al. 2017, MNRAS, 465, 4795 [NASA ADS] [CrossRef] [Google Scholar]
 Bouret, J.C., Hillier, D. J., Lanz, T., & Fullerton, A. W. 2012, A&A, 544, A67 [Google Scholar]
 Carciofi, A. C., Domiciano de Souza, A., Magalhães, A. M., Bjorkman, J. E., & Vakili, F. 2008, ApJ, 676, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Chaboyer, B., & Zahn, J.P. 1992, A&A, 253, 173 [Google Scholar]
 Che, X., Monnier, J. D., Zhao, M., et al. 2011, ApJ, 732, 68 [NASA ADS] [CrossRef] [Google Scholar]
 Chieffi, A., & Limongi, M. 2013, ApJ, 764, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Cohen, D., Sundqvist, J., & Leutenegger, M. 2013, in Massive Stars: From Alpha to Omega, 36 [Google Scholar]
 Crowther, P. A., Lennon, D. J., & Walborn, N. R. 2006, A&A, 446, 279 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Curé, M. 2004, ApJ, 614, 929 [NASA ADS] [CrossRef] [Google Scholar]
 Curé, M., Rial, D. F., & Cidale, L. 2005, A&A, 437, 929 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Decressin, T., Mathis, S., Palacios, A., et al. 2009, A&A, 495, 271 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Domiciano de Souza, A., Kervella, P., Moser Faes, D., et al. 2014, A&A, 569, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eggenberger, P., Meynet, G., Maeder, A., et al. 2008, Ap&SS, 316, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Espinosa Lara, F., & Rieutord, M. 2011, A&A, 533, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Espinosa Lara, F., & Rieutord, M. 2013, A&A, 552, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Friend, D. B., & Abbott, D. C. 1986, ApJ, 311, 701 [NASA ADS] [CrossRef] [Google Scholar]
 Gagnier, D., Rieutord, M., Charbonnel, C., Putigny, B., & Espinosa Lara, F. 2019, A&A, 625, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gayley, K. G. 1995, ApJ, 454, 410 [NASA ADS] [CrossRef] [Google Scholar]
 Georgy, C., Meynet, G., & Maeder, A. 2011, A&A, 527, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Georgy, C., Ekström, S., Granada, A., et al. 2013, A&A, 553, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Glatzel, W. 1998, A&A, 339, L5 [NASA ADS] [Google Scholar]
 Granada, A., & Haemmerlé, L. 2014, A&A, 570, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Granada, A., Ekström, S., Georgy, C., et al. 2013, A&A, 553, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Grevesse, N., & Noels, A. 1993, in Origin and Evolution of the Elements, eds. N. Prantzos, E. VangioniFlam, & M. Casse, 15 [Google Scholar]
 Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368 [NASA ADS] [CrossRef] [Google Scholar]
 Hervé, A., Rauw, G., & Nazé, Y. 2013, A&A, 551, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jackson, S., MacGregor, K. B., & Skumanich, A. 2005, ApJS, 156, 245 [NASA ADS] [CrossRef] [Google Scholar]
 Krtička, J., Owocki, S. P., & Meynet, G. 2011, A&A, 527, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kudritzki, R.P., & Puls, J. 2000, ARA&A, 38, 613 [NASA ADS] [CrossRef] [Google Scholar]
 Kudritzki, R. P., Pauldrach, A., & Puls, J. 1987, A&A, 173, 293 [NASA ADS] [Google Scholar]
 Kudritzki, R.P., Lennon, D. J., & Puls, J. 1995, in Science with the VLT, ed. J. R. Walsh, & I. J. Danziger (Berlin: SpringerVerlag), 246 [CrossRef] [Google Scholar]
 Lamers, H. J. G., & Pauldrach, A. W. A. 1991, A&A, 244, L5 [NASA ADS] [Google Scholar]
 Lamers, H. J. G. L. M., Snow, T. P., & Lindholm, D. M. 1995, ApJ, 455, 269 [NASA ADS] [CrossRef] [Google Scholar]
 Langer, N. 1997, in Luminous Blue Variables: Massive Stars in Transition, eds. A. Nota, & H. Lamers, ASP Conf. Ser., 120, 83 [Google Scholar]
 Langer, N. 1998, A&A, 329, 551 [NASA ADS] [Google Scholar]
 Lau, H. H. B., Potter, A. T., & Tout, C. A. 2011, MNRAS, 415, 959 [NASA ADS] [CrossRef] [Google Scholar]
 Leutenegger, M. A., Cohen, D. H., Sundqvist, J. O., & Owocki, S. P. 2013, ApJ, 770, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Lignieres, F., Catala, C., & Mangeney, A. 2000, ArXiv eprints [arXiv:astroph/0002026] [Google Scholar]
 Maeder, A. 1999, A&A, 347, 185 [NASA ADS] [Google Scholar]
 Maeder, A. 2009, Physics, Formation and Evolution of Rotating stars (Berlin: Springer) [Google Scholar]
 Maeder, A., & Meynet, G. 2000, A&A, 361, 159 [NASA ADS] [Google Scholar]
 Maeder, A., & Meynet, G. 2010, New Astron. Rev., 54, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Maeder, A., & Meynet, G. 2015, in New Windows on Massive Stars, eds. G. Meynet, C. Georgy, J. Groh, & P. Stee, IAU Symp., 307, 9 [NASA ADS] [Google Scholar]
 Maeder, A., & Zahn, J. P. 1998, A&A, 334, 1000 [Google Scholar]
 Markova, N., & Puls, J. 2008, A&A, 478, 823 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marques, J. P., Goupil, M. J., Lebreton, Y., et al. 2013, A&A, 549, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Meynet, G., & Maeder, A. 1997, A&A, 321, 465 [NASA ADS] [Google Scholar]
 Meynet, G., & Maeder, A. 2000, A&A, 361, 101 [NASA ADS] [Google Scholar]
 Meynet, G., Chomienne, V., Ekström, S., et al. 2015, A&A, 575, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Monnier, J. D., Zhao, M., Pedretti, E., et al. 2007, Science, 317, 342 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Müller, P. E., & Vink, J. S. 2014, A&A, 564, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Najarro, F., Hanson, M. M., & Puls, J. 2011, A&A, 535, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Owocki, S. P., & Gayley, K. G. 1997, in Luminous Blue Variables: Massive Stars in Transition, eds. A. Nota, & H. Lamers, ASP Conf. Ser., 120, 121 [NASA ADS] [Google Scholar]
 Owocki, S. P., Cranmer, S. R., & Gayley, K. G. 1996, ApJ, 472, L115 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P., Gayley, K. G., & Cranmer, S. R. 1998a, in Properties of Hot Luminous Stars, ed. I. Howarth, ASP Conf. Ser., 131, 237 [Google Scholar]
 Owocki, S. P., Cranmer, S. R., & Gayley, K. G. 1998b, in B[e] stars, eds. A. M. Hubert, & C. Jaschek, Astrophys. Space Sci. Lib., 233, 205 [NASA ADS] [CrossRef] [Google Scholar]
 Pauldrach, A., Puls, J., & Kudritzki, R. P. 1986, A&A, 164, 86 [NASA ADS] [Google Scholar]
 Pauldrach, A. W. A., Kudritzki, R. P., Puls, J., & Butler, K. 1990, A&A, 228, 125 [NASA ADS] [Google Scholar]
 Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Pelupessy, I., Lamers, H. J. G. L. M., & Vink, J. S. 2000, A&A, 359, 695 [NASA ADS] [Google Scholar]
 Petrenz, P., & Puls, J. 1996, A&A, 312, 195 [NASA ADS] [Google Scholar]
 Petrenz, P., & Puls, J. 2000, A&A, 358, 956 [NASA ADS] [Google Scholar]
 Petrov, B., Vink, J. S., & Gräfener, G. 2016, MNRAS, 458, 1999 [NASA ADS] [CrossRef] [Google Scholar]
 Porter, J. M., & Rivinius, T. 2003, PASP, 115, 1153 [NASA ADS] [CrossRef] [Google Scholar]
 Przybilla, N., Nieva, M.F., & Butler, K. 2008, ApJ, 688, L103 [NASA ADS] [CrossRef] [Google Scholar]
 Puls, J., Kudritzki, R.P., Herrero, A., et al. 1996, A&A, 305, 171 [NASA ADS] [Google Scholar]
 Puls, J., Springmann, U., & Lennon, M. 2000, A&AS, 141, 23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Puls, J., Sundqvist, J. O., & Markova, N. 2015, in New Windows on Massive Stars, eds. G. Meynet, C. Georgy, J. Groh, & P. Stee, IAU Symp., 307, 25 [Google Scholar]
 Puls, J., Vink, J. S., & Najarro, F. 2008, A&ARv, 16, 209 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rauw, G., Hervé, A., Nazé, Y., et al. 2015, A&A, 580, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Renzo, M., Ott, C. D., Shore, S. N., & de Mink, S. E. 2017, A&A, 603, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rieutord, M. 2016, in Cartography of the Sun and the Stars, eds. J. P. Rozelot, & C. Neiner (Berlin Springer Verlag), Lect. Notes Phys., 914, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Rieutord, M., & Beth, A. 2014, A&A, 570, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rieutord, M., Espinosa Lara, F., & Putigny, B. 2016, J. Comput. Phys., 318, 277 [NASA ADS] [CrossRef] [Google Scholar]
 Rivinius, T., Carciofi, A. C., & Martayan, C. 2013, A&ARv, 21, 69 [Google Scholar]
 Roxburgh, I. W. 2004, A&A, 428, 171 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Runacres, M., & Blomme, R. 1994, Ap&SS, 216, 69 [NASA ADS] [CrossRef] [Google Scholar]
 SimónDíaz, S., & Herrero, A. 2007, A&A, 468, 1063 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Smith, N. 2014, ARA&A, 52, 487 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Sundqvist, J. O., & Owocki, S. P. 2013, MNRAS, 428, 1837 [NASA ADS] [CrossRef] [Google Scholar]
 Sundqvist, J. O., Puls, J., Feldmeier, A., & Owocki, S. P. 2011, A&A, 528, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 1999, A&A, 350, 181 [Google Scholar]
 Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vink, J. S., Brott, I., Gräfener, G., et al. 2010, A&A, 512, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 von Zeipel, H. 1924, MNRAS, 84, 665 [NASA ADS] [CrossRef] [Google Scholar]
 Yoon, S.C., & Langer, N. 2005, A&A, 443, 643 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zahn, J.P. 1992, A&A, 265, 115 [NASA ADS] [Google Scholar]
 Zhao, M., Monnier, J. D., Pedretti, E., et al. 2009, ApJ, 701, 209 [NASA ADS] [CrossRef] [Google Scholar]
 Zickgraf, F. J., Wolf, B., Stahl, O., Leitherer, C., & Appenzeller, I. 1986, A&A, 163, 119 [NASA ADS] [Google Scholar]
 Zickgraf, F.J., Wolf, B., Stahl, O., & Humphreys, R. M. 1989, A&A, 220, 206 [Google Scholar]
 Zorec, J., Rieutord, M., Espinosa Lara, F., et al. 2017, A&A, 606, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
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,
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
where τ_{s} ≃ 2/3 is the Rosseland mean optical depth at the photosphere. Assuming the power law and the ideal gas equation of state, the previous expression leads to
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 largescale 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,
where ϕ is is the gravitational potential; the continuity equation,
the momentum equation,
where F_{visc} is the viscous force; and the heat balance equation,
and
This last equation assumes an efficient convection in convective cores, as can be shown with the mixinglength model.
These equations are completed by boundary conditions that require that (i) the gravitational potential ϕ vanishes at infinity, (ii) velocity fields meet stressfree 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 microphysics, 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
Summary of the main results for ESTER 2D models of a 15 M_{⊙} star with Z = 0.02 at ZAMS.
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 g_{eff}(π/2)=0. 

In the text 
Fig. 2. Relative difference between the radiative flux of the ωmodel and the ESTER model as a function of colatitude for a 15 M_{⊙} (top) and a 40 M_{⊙} (bottom) ZAMSstar 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 X_{core}/X_{0} = 0.5. X_{0} is the initial hydrogen mass fraction at ZAMS. The minimum of each curve corresponds to a sign change of F_{ESTER} − F_{ω}. 

In the text 
Fig. 3. Rotationdependent 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 colatitude 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 T_{eff} 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 T_{eff} 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 (bottom) as a function of colatitude θ for a 15 M_{⊙} ESTER 2Dmodel at ZAMS with Z = 0.02 and various angular velocity ratios ω. The main bistability limit is reached near the equator for ω ≳ 0.85, and induces a strong massflux 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 massloss rate Ṁ (in M_{⊙} yr^{−1}) for a 15 M_{⊙} (red) and a 10 M_{⊙} (green) ESTER model against the corresponding mean effective temperature , at ZAMS with Z = 0.02 and for ω ∈ [0, 1[. Arrows indicate the direction of increasing ω. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.