Critical angular velocity and anisotropic mass loss of rotating stars with radiation-driven winds

The understanding of the evolution of early-type stars is tightly related to that of the effects of rapid rotation. For massive stars, rapid rotation combines with their strong radiation-driven wind. 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? To investigate fast rotation, which makes stars oblate, we used the 2D ESTER models and a simplified approach, the $\omega$-model, which gives the latitudinal dependence of the radiative flux in a centrifugally flattened radiative envelope. We find that radiative acceleration only mildly influences the critical angular velocity, at least for stars with masses lower than 40 Msun. 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 so-called bi-stability jump) included in the $\dot{M}-T_{\rm eff}$ relation of 1D models means that the mass flux of a fast-rotating star is controlled by either a single wind or a two-wind regime. Mass and angular momentum losses are strong around the equator if the star is in the two-wind regime. We also show that the difficulty of selecting massive stars that are viewed pole-on makes detecting the discontinuity in the relation between mass loss and effective temperature also quite challenging.


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. , 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. 2016FRANEC, Chieffi & Limongi 2013; CES-TAM, 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 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 con-sidered 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 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;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 & Beth2014).
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).

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 flux-weighted opacity 1 , 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  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 break-up limit (or the Ω-limit) when radiative acceleration can be neglected in total gravity. In this equation, G is the 1 Strictly speaking, κ is the mass absorption coefficient, κ = µ/ρ, with opacity µ = λ −1 , where λ is the mean free path of photons, and density ρ. However, in the following, we use, as is frequently done in the current context, the term "opacity". 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 (1997Langer ( , 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  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 where χ = 4acT 3 /(3κρ) is the radiative conductivity. Additionally, assuming solid-body 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 solid-body 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.
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 . These limitations prompt us to revisit this limit with ESTER 2D models.

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 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 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 cos ψ + ln tan(ψ/2) = 1 3 ω 2 r 3 cos 3 θ + cos θ + ln tan(θ/2). (14) 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

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 rotation-dependent 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.

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.  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).

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 g tot (π/2) = g eff (π/2) + g rad (π/2), where g rad (π/2) = − κ(π/2)L 4πcGM 1 − ω 2 −2/3 g eff (π/2), 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 Ω 2 k R eq as g tot (π/2) = ω 2 + Γ eq (1 − ω 2 ) 1/3 − 1, 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. 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 ω.

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 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. 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 4 eff /g eff . 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 .  came to the same conclusion, but with a different expression for critical angular velocity, namely Ω c ∝ Ω k 1 − Γ 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  for critical angular velocity in this regime.

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.
A88, page 4 of 14 2. This balance of forces is always first reached at the equator when Ω/Ω k increases because the ratio T 4 eff /g 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.

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 eff ∝ g β eff . 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 4 eff /g eff at the equator when criticality is approached.

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 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 as a function of co-latitude. The result is shown in Fig. 2. 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.

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, 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 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).
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 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 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 .
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 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 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.
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.

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 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 1979or 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).

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, and the radial momentum equation, 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, where g line rad is the line-driven 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 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, 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 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) 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 . 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 line-driving 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 mass-loss rate of a non-rotating star, namelẏ The global mass loss (without rotation) thus scales aṡ which is the basis of the wind momentum-luminosity relation (Kudritzki et al. 1995;Puls et al. 1996).

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 µ * = √ 1 − u 2 . 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 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 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 mass-loss rate. Thus we usė 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 11 cm −3 and W ≡ 0.5(1 − √ 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 α ≡ α − δ . Finally, we obtain the modified local mass-flux in the non-rotating 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.

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 .
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. 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 jump eff 20 kK, while Vink et al. (2001) predicted anṀ-jump by a factor ∼10 located at T jump eff 25 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 eff > T jump eff 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;Puls et al. 2015) and/or to the outdated solar mixture used in 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 B T 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, 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.

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 readṡ 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 3 (e.g. Owocki & Gayley 1997;Owocki et al. 1998a;Petrenz & Puls 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 ω.

Latitudinal variations in mass and angular momentum loss
We computed the local mass-fluxṁ(θ) as well as the local angular momentum fluẋ with ESTER 2D models and prescription (48). Figure 8 shows the local mass-fluxṁ 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 T eff (θ jump ) = T jump eff 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 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 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. (1986Zickgraf et al. ( , 1989, Lamers & Pauldrach (1991), Owocki & Gayley (1997), Owocki et al. (1998b), and Pelupessy et al. (2000).
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 ω.

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/∂θ . 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 T L = L/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).
With increasing ω, the total angular momentum of the star L increases, and even though the global mass-loss rateṀ decreases in the SWR, the global loss of angular momentumL also increases in this regime. This is simply becauseL increases for increasing ω. It is even more interesting, however, that the timescale of angular momentum loss T L = L/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 L rapidly decreases as ω approaches unity. Both the strong increase inṀ and decrease in T L can be explained by the increasing stellar surface fraction where the effective temperature is lower than T jump eff 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 mass-loss rate if the mean effective temperature of the star is lower than T jump eff ∼ 22.5-25 kK, according to Vink et al. (2001). In the present work however, 2D models reach the bi-stability limit Notes. The first column reports the ratio of equatorial angular velocity to Keplerian angular velocity ω, the second column reports the critical angular velocity ratio Ω eq /Ω c , the third column is the global mass-loss rate logṀ, and the last column gives the ratio between total angular momentum and angular momentum loss rate log T L = log L/L.Ṁ is in M yr −1 and T L in yr. (a) Star in the TWR.
if the local effective temperature somewhere on the stellar surface is lower than T jump eff . This difference has two consequences. Firstly, 2D models can reach the bi-stability limit even with an average effective temperatures higher than T jump eff . 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 surfaceaveraged effective temperature T eff of the model. In these models, T eff is greater than T jump eff 22.8 kK for all ω. Thus, equivalent 1D models would just have ignored the bi-stability jump.
Secondly, in 2D models the surface fraction where the effective temperature is lower than T jump eff monotonically increases with increasing ω. This results in a gradual variation inṀ (and T L ) with T 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 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 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 v sin i to select either slowly rotating stars or stars that are viewed pole-on. Obviously, the precise determination of v sin i for hot massive stars is a challenge in particular because rotational effects are mixed with other linebroadening effects such as the so-called macro-turbulence (e.g. Simón-Díaz & Herrero 2007).

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. 2000Puls et al. , 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.

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 . 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 . 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 (1997Langer ( , 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 followup 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.