© The Authors 2023

1. Introduction

All massive stars suffer from the effects of stellar winds. For O-type stars, the winds can be so strong that a significant portion of the star evaporates and the evolutionary pathway is altered dramatically, for example forming a Wolf-Rayet star (Maeder & Meynet 1987; de Koter et al. 1997). In contrast, lower-mass stars typically lose only a negligible fraction of their mass to winds. However even a non-magnetic wind carries angular momentum away from a star, and a star’s spin evolution can even be affected by weak winds (Langer 1998).

In the simplest sense, the density and velocity structure and thus also the mass-loss rate of a radiation-driven wind is determined by the opposing effects of gravity and radiative acceleration. Gravity serves to bind material to the stellar surface, while radiation, through both continuum and line opacities, provides a force to overcome gravity (Castor et al. 1975; Pauldrach et al. 1986). Rotation directly affects both the gravitational field strength and the radiation field (von Zeipel 1924), with both varying over the stellar surface, in turn resulting in an anisotropic wind (Poe & Friend 1986; Cranmer & Owocki 1995; Curé & Rial 2004).

For an anisotropic wind, attention needs to be paid to angular momentum loss since mass lost at the equator carries a larger specific angular momentum than mass lost near the poles, especially so for stars that are significantly deformed from sphericity. Owing to internal structural changes, stars born with moderate rotation may evolve to become extremely fast rotators (Hastings et al. 2020b), so the effects of rotation on stellar winds have the potential to be relevant to a large portion of stars.

For a number of decades, massive star modelling efforts (Heger et al. 2000; Brott et al. 2011; Paxton et al. 2013) have described the effects of rotation on mass loss by increasing the mass-loss rate of an equivalent non-rotating star by a factor that depends on the rotation rate (Friend & Abbott 1986). Such a formulation is lacking due to two issues. Firstly, it is assumed that, independent of the wind recipe used, rotation always increases mass-loss rates by the same relative amount. This is not a fair assumption because two of the major effects of rotation on the surface of a star are to weaken the gravitational field and to reduce the surface-averaged effective temperature (von Zeipel 1924). These two effects generally, though not always, serve to counteract each other, with winds being enhanced by weaker gravities but diminished by lower effective temperatures. It is unclear which effect dominates. Both the dependence of the mass-loss rate on temperature and the assumed temperature profile across the stellar surface (gravity darkening law) will govern whether rotation enhances or reduces mass loss, meaning that the enhancement ought to be model dependant (cf. Müller & Vink 2014).

Secondly, some of the mathematical functions used to provide the mass-loss enhancement diverge as the star approaches the critical velocity. While this behaviour is used in stellar evolution calculations merely to prevent models from exceeding critical rotation (Heger et al. 2000; Petrovic et al. 2005), it is unphysical not least because it is usually¹ only material at the equator which achieves the critical rotation velocity, and strictly the equator covers an infinitesimally small surface, while gravity does manage to keep the star bound over the rest of the surface.

Angular momentum loss from massive stars plays a role in several active research topics such as the study of Be stars (Curé 2004; Curé et al. 2005; Ekström et al. 2008; Hastings et al. 2020b); the occurrence of chemically homogeneous evolution, relevant to double black-hole mergers (Marchant et al. 2016) and gamma-ray burst progenitors (Yoon et al. 2006; Aguilera-Dena et al. 2020); wind-driven orbital evolution in massive binary stars (MacLeod & Loeb 2020; Sen et al. 2022); and of course the rotation rates of stars in general. Improved modelling of the winds of rotating stars would be beneficial to the advancement our understanding of stellar physics.

Various studies concerning winds from rotating massive stars have been performed (e.g., Poe & Friend 1986; Poe et al. 1989; Owocki et al. 1994; Cranmer & Owocki 1995; Pelupessy et al. 2000; Petrenz & Puls 2000; Curé et al. 2012; Müller & Vink 2014; Gagnier et al. 2019b), although the results of which have not been adopted for use in stellar evolution codes. Our aim is to provide a simple and easily implementable scheme that improves upon the popular rotationally enhanced mass-loss schemes. We applied our scheme to the one-dimensional stellar evolution code MESA and provide the files necessary to compute models using it².

The structure of this paper is as follows. Section 2 details the derivation of our stellar wind prescription. In Sect. 3 we compare the results of our new wind model to the commonly used rotationally enhanced mass-loss prescription. A brief discussion of uncertainties is given in Sect. 4. Section 6 hosts a comparison of our results to more sophisticated approaches. Lastly, our conclusions are put forward in Sect. 7.

2. Method

2.1. Anisotropic wind model

Our basic philosophy is to apply a one-dimensional wind recipe to every point on the surface of a rotating star. For every point on the stellar surface, the given wind recipe uses the local physical conditions to provide a surface mass-flux, which when integrated results in global mass and angular momentum loss rates. We shall now determine the surface properties of a rotating star.

2.1.1. Surface properties of a rotating star

A rotating star with mass M, polar radius R_p, equatorial radius R_e, luminosity L is assumed to be rotating rigidly with angular velocity Ω. In reasonable agreement with detailed stellar models (Maeder 2009), we assume that the polar radius is not affected by rotation. The contribution of radiative acceleration to the total gravity shall be ignored, as we focus primarily on stars with luminosities below the Eddington luminosity. The critical velocity, or break-up velocity is then the Keplerian angular velocity, at which the gravitational force matches the centrifugal force at the equator and reads

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\text{Kep}}=\sqrt{\frac{\mathit{GM}}{R_{\mathrm{e}}^3}}\end{array}$$(1)

and the fraction of Keplerian angular velocity is denoted as

$$\begin{array}{c}\hfill \omega =\frac{\mathrm{\Omega }}{{\mathrm{\Omega }}_{\text{Kep}}}.\end{array}$$(2)

Recent two-dimensional models of rotating stars suggest that the rotation velocity at which material becomes unbound from the stellar surface is very close to the Keplerian velocity (Gagnier et al. 2019a). However, these models only cover two separate values of the stellar mass (15 and 40 M_⊙) at one point in their evolution. Therefore we cannot exclude that for very luminous stars, radiation might play a significant role in unbinding material from the surface and thus reducing the critical rotation velocity. This issue is discussed further in Sect. 2.2.

We note that several different working definitions of the critical velocity velocity exist (see discussion in Sect. 2.3.1 of Rivinius et al. 2013 and their Eqs. (3) and (4)). Our choice is made to be consistent with the stellar evolution code MESA (Paxton et al. 2019).

In the co-rotating frame of a rotating star, the centrifugal force is perpendicular to the rotation axis, which causes the effective surface gravity to to become latitude-dependant. Following from the varying surface gravity, effective temperature also varies across the surface (von Zeipel 1924). Also effected is the star’s shape, evidenced by a bulging equator. These three effects shall now be quantified in order.

As massive stars are centrally condensed, the use of the Roche potential is justified (Collins 1965; Rieutord 2016). We define the effective surface gravity as the sum of self-gravitation and centrifugal forces which is

$$\begin{array}{c}\hfill {\mathit{g}}_{\text{eff}}(\theta )=(-\frac{\mathit{GM}}{r{(\theta )}^2}\text{sin}(\theta )+{\mathrm{\Omega }}^{2}r(\theta )\text{sin}(\theta ))\mathit{x}-\frac{\mathit{GM}}{r{(\theta )}^2}\text{cos}(\theta )\mathit{z},\end{array}$$(3)

where x and z are the Cartesian unit vectors, perpendicular and parallel to the rotation axis respectively. The radial co-ordinate is designated r and θ the co-latitude. The magnitude of the surface gravity is then found to be

$$\begin{array}{cc}\hfill |{\mathit{g}}_{\text{eff}}(\theta )|& =\left(\frac{\mathit{GM}}{R_{\mathrm{p}}^2}\right){\delta }^{-2}[{\delta }^4{\left(\frac{r(\theta )}{R_{\mathrm{p}}}\right)}^{-4}\hfill \\ \hfill & +{\omega }^4{\delta }^{-2}{\left(\frac{r(\theta )}{R_{\mathrm{p}}}\right)}^2{\text{sin}}^2(\theta )-2{\omega }^2{\text{sin}}^2(\theta )\delta {\left(\frac{r(\theta )}{R_{\mathrm{p}}}\right)}^{-1}{]}^{\frac{1}{2}},\hfill \end{array}$$(4)

where δ is the ratio of equatorial and polar radii, $\delta =\frac{R_{\mathrm{e}}}{R_{\mathrm{p}}}$. Because we assume that every point on the surface can be treated as an equivalent non-rotating star (i.e., we wish to reduce the two-dimensional problem of a rotating star to one dimension), and that the flux vector is nearly perfectly aligned with the gravity vector in a rotating star (Espinosa Lara & Rieutord 2011), the magnitude of the gravity vector is the quantity of interest, not the gravitational field strength in the radial direction.

The local effective temperature is defined using the local flux, F(θ), and the Stefan-Boltzmann constant, σ as

$$\begin{array}{c}\hfill T_{\mathrm{eff}}{(\theta )}^4=|\mathit{F}(\theta )|/\sigma .\end{array}$$(5)

The effective temperature profile is given by the model of Espinosa Lara & Rieutord (2011), which assumes a Roche potential and that the flux at the surface of a star is well approximated by

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

which requires the energy flux to be anti-parallel to the effective gravity. This condition is fulfilled in stars with convective envelopes and is also valid to within a very fine tolerance in stars with radiative envelopes (Espinosa Lara & Rieutord 2011). The function f(r, θ) is found by demanding that no heat is generated in the envelope (i.e., ∇F = 0) and reads

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

where ϑ is the solution to

$$\begin{array}{c}\hfill \text{cos}\vartheta +\text{ln}\text{tan}\frac{\vartheta }{2}=\frac{1}{3}{\omega }^2{\left(\frac{r}{R_{\mathrm{e}}}\right)}^3{\text{cos}}^3\theta +\text{cos}\theta +\text{ln}\text{tan}\frac{\theta }{2}.\end{array}$$(8)

Alternative gravity darkening laws are available (Slettebak 1949; Lucy 1967; Lovekin et al. 2006; Lipatov & Brandt 2020). We note that the gravity darkening model of Espinosa Lara & Rieutord (2011) predicts that the equatorial flux of a critically rotating star is zero, which might be unphysical.

Lastly, the radial profile can be determined from the Roche equipotential surface (Appendix A) to be

$$\begin{array}{c}\hfill \frac{r(\omega ,\theta )}{R_{\mathrm{p}}}=(2+{\omega }^2)\sqrt{\frac{2+{\omega }^2}{3{\omega }^2{\text{sin}}^2(\theta )}}\text{cos}[\frac{1}{3}\text{arccos}\left(\frac{3}{2+{\omega }^2}\sqrt{\frac{3{\omega }^2{\sin }^2(\theta )}{2+{\omega }^2}}\right)+\frac{\pi }{3}].\end{array}$$(9)

2.1.2. Mass and angular momentum flux

To quantify the wind over the stellar surface, we shall use the mass-loss rate per unit surface area, or mass-flux, ṁ(θ), which is related to the total mass-loss rate, $\dot{\mathcal{M}}$ via

$$\begin{array}{c}\hfill \dot{\mathcal{M}}={\displaystyle \int }\dot{m}(\theta )\mathrm{d}S,\end{array}$$(10)

where dS represents the infinitesimal surface element. Knowledge of the star’s shape allows us to compute the above integral as (cf. Gagnier et al. 2019b)

$$\begin{array}{c}\hfill \dot{\mathcal{M}}=2\pi {\displaystyle \int }\dot{m}(\theta )r^2(\theta )\sqrt{1+\frac{1}{r^2(\theta )}{\left(\frac{\partial r}{\partial \theta }\right)}^2}\text{sin}\theta \mathrm{d}\theta ,\end{array}$$(11)

where r(θ) is given by Eq. (9).

The local angular momentum flux is defined as

$$\begin{array}{c}\hfill \dot{\ell }(\theta )=\dot{m}\mathrm{\Omega }{r}^2(\theta ){\text{sin}}^2\theta \end{array}$$(12)

and the global angular momentum loss rate is found by integrating again over the stellar surface as

$$\begin{array}{c}\hfill \dot{\mathcal{L}}=2\pi {\displaystyle \int }\dot{\ell }(\theta )r^2(\theta )\sqrt{1+\frac{1}{r^2(\theta )}{\left(\frac{\partial r}{\partial \theta }\right)}^2}\text{sin}\theta \mathrm{d}\theta .\end{array}$$(13)

2.1.3. Determining surface mass flux

Calculating the surface mass flux of a rotating massive star requires not only knowledge of the general mechanics of radiation-driven winds but also of several rotation specific phenomena and their interplay in driving a wind. As of yet, general mass-loss recipes exist only for non-rotating stars, and even those differ significantly depending on methods and assumptions. It is felt that although the use of a non-rotating wind recipe cannot capture the fine details of physical processes in rotating stars, their use in describing rotating star winds is still beneficial and above all represents an improvement over the almost exclusively used rotationally enhanced mass-loss scheme.

Using the effective surface gravity profile, effective temperature profile and surface shape of a rotating star, we may assign an equivalent non-rotating star to each co-latitude of the rotating star, for which the mass-loss rate is given by a chosen recipe. This equivalent non-rotating star is defined to have the same radius, effective temperature and surface gravity as a given latitude on the rotating star. The surface mass-flux is then, in general

$$\begin{array}{c}\hfill \dot{m}(\theta )=\frac{\dot{M}(|{\mathit{g}}_{\text{eff}}(\theta )|,T_{\mathrm{eff}}(\theta ),r(\theta ),\dots )}{4\pi r{(\theta )}^2},\end{array}$$(14)

where Ṁ is the function provided by the non-rotating wind recipe. The only requirement for the recipe is that it is a function of, or can be manipulated to be a function of, at least the effective surface gravity, effective temperature and radius.

For the calculations in this work, we shall use the mass-loss recipe of Vink et al. (2001), although in principle any recipe can be used. Here the mass-loss rate is a function of the stellar mass, effective temperature, luminosity and metallicity, Z. For our purposes, we first need to modify the input parameters of the recipe.

The mass and luminosity of a non-rotating star can be described using the effective surface gravity and effective temperature, provided the radius is known. This means that at each latitude of a rotating star, an equivalent non-rotating star would have a different mass (following from the radius and surface gravity of the rotating star) and a different luminosity (following from the effective temperature and radius). To account for this, the mass-dependence must be expressed instead in terms of the surface gravity and radius. As luminosity is determined by the Stefan-Boltzmann law, an effective luminosity for each colatitude on a rotating star can be defined as

$$\begin{array}{c}\hfill L_{\text{eff}}(\theta )=4\pi \sigma r{(\theta )}^2{T}_{\mathrm{eff}}{(\theta )}^4.\end{array}$$(15)

This equivalent luminosity represents the luminosity that a non-rotating star, with equal surface properties of a given colatitude, would have. It is therefore this quantity that must be used in place of the luminosity in the mass-loss recipe, which becomes

$$\begin{array}{c}\hfill \dot{m}(\theta )=\frac{\dot{M}(\frac{|{\mathit{g}}_{\text{eff}}{(\theta )|r(\theta )}^2}{G},T_{\mathrm{eff}}(\theta ),L_{\text{eff}}(\theta ),Z)}{4\pi r{(\theta )}^2}.\end{array}$$(16)

2.2. Critical rotation velocity

For a star, there exists a critical (or break-up) rotation velocity at which material becomes unbound from the stellar surface. Although a simple concept, there are several nuances which shall be discussed here. In this work we assume that the critical velocity is attained when the centrifugal and gravitational forces balance, however in general this is likely only an approximation.

In massive stars the force from radiation itself contributes to the force balance, and thus has been proposed to reduce the critical velocity (Langer 1997). The acceleration produced by radiation is proportional to the flux and opacity, which are both effected by rotation. As first argued by Glatzel (1998), when a luminous star rotates very quickly, gravity darkening causes the equatorial flux to weaken strongly, suggesting that the appropriate limit is the Keplerian one. Although analysis by Maeder & Meynet (2000) determined that below a threshold luminosity (around 60% of the Eddington luminosity), the radiation force indeed plays no role in unbinding material from the surface, the issue is still not clear cut, as discussed in the following.

Gravity darkening is traditionally described by Von Zeipel’s Law, which states that effective temperature is proportional to effective gravity to the power of β, with β = 0.25. Interferometric observations of rapidly rotating stars have demonstrated that gravity darkening is not as strong as predicted by Von Zeipel’s Law, with lower β values for faster rotators (Monnier et al. 2007; Zhao et al. 2009; Che et al. 2011; Domiciano de Souza et al. 2014). These observations are generally consistent with two-dimensional numerical models (Espinosa Lara & Rieutord 2013) and analytic gravity darkening models (Espinosa Lara & Rieutord 2011), but one star, β Cassiopeiae, appears to exhibit much weaker gravity darkening than expected (Che et al. 2011), perhaps exposing weaknesses in our understanding of gravity darkening. Weaker gravity darkening would result in a stronger radiative force at the equator, hence helping to reduce the critical velocity.

The surface opacity of a rotating star is also uncertain. Maeder & Meynet (2000) assumed that the region with the highest opacity would be the equator, as this is the coldest part of the surface. However, the centrifugal force also causes a decrease in the matter density at the equator, consequently decreasing the opacity. Two-dimensional numerical models of stars on the zero-age-main-sequence suggest that the effect of decreasing density dominates, thus a fast rotating star is predicted to have a lower equatorial opacity than an equivalent non-rotating star (Gagnier et al. 2019a), meaning that radiative acceleration is unable to contribute to the force balance. However there may be some exceptions. Firstly stars may suffer the effects of opacity bumps caused by recombination of certain species (notably hydrogen, helium and iron; Iglesias & Rogers 1996), that could drastically alter the opacity profile over the surface of the star. Secondly in very hot stars, where the opacity is dominated by electron scattering, the surface opacity is largely independent of temperature and thus spatially uniform. Such a case would need careful study to determine whether the break-up velocity is affected.

Classical Be stars are fast rotators with a decretion disc, which are believed to be typically rotating at approximately 70% of the Keplerian velocity (Porter 1996; Rivinius et al. 2013; Zorec et al. 2016; Dufton et al. 2022), and in some cases even lower (Huang et al. 2010; Zorec et al. 2016). It may be reasonably supposed that an outflowing disc will affect the structure of its host star, such that the surface flux and opacities may be different when a disc is present, thus altering the break-up velocity. There is evidence to suggest that the threshold rotation rate for the Be phenomenon, assumed to be the true break-up velocity, varies with effective temperature (Cranmer 2005; Huang et al. 2010), with hotter Be stars rotating more slowly than their cool counterparts. It is well understood that hotter stars are closer to the Eddington limit, which may suggest that indeed in the hotter Be stars, radiative acceleration does play a role in unbinding material.

2.3. Numerical method

In order to investigate the effect of our prescription on the evolution of both mass and angular momentum loss rates of rotating stars, we employ the one-dimensional detailed stellar evolution code MESA (Paxton et al. 2019), version 12115. The files required to compute models presented in this work are available online³. The adopted physics is largely identical, except for the stellar winds, to that of Brott et al. (2011) and implemented in MESA as by Schootemeijer et al. (2019). The models include internal angular momentum transport achieved by magnetic torques (Spruit 2002) which enforce near solid-body rotation during most of the main-sequence evolution.

We run two sets of models, one using the rotationally enhanced mass-loss scheme as it is usually implemented in MESA (named the standard scheme), where the mass-loss rates are first calculated using the recipe of Vink et al. (2001) and then following Friend & Abbott (1986)⁴ increased by a factor of

$$\begin{array}{c}\hfill {(1-\mathrm{\Omega }\sqrt{\frac{R_{\mathrm{e}}^3}{GM(1-\mathrm{\Gamma })}})}^{-0.43}.\end{array}$$(17)

The second set uses mass-loss rates set by the method outlined in Sect. 2.1 and is named the local scheme. Both sets rely on the wind mass-loss recipe of Vink et al. (2001). This wind recipe includes the bi-stability jump effect (first found by Pauldrach & Puls 1990), where mass-loss rates are theorised to increase dramatically during the transition to temperatures cooler than approximately 22 kK owing to the recombination of Fe IV in the atmosphere (Vink et al. 1999). The impact of the bi-stability jump on mass-loss rates is not certain, with Björklund et al. (2023) noting that ‘the drastic Ṁ increase found in earlier models in this region might simply be an artefact of not being dynamically consistent around the sonic point, and not allowing properly for the feedback between radiative and velocity acceleration’. The quantitative behaviour of models near the jump is also contested (Markova & Puls 2008; Vink 2018; Krtička et al. 2021). We stress that our method is not confined to a particular mass-loss recipe and that several others could be used, for example those of Kudritzki et al. (1989), Sundqvist et al. (2019), Björklund et al. (2023).

For the standard scheme, stellar winds are assumed to be isotropic with the angular momentum loss $\dot{\mathcal{L}}$, given by

$$\begin{array}{c}\hfill \dot{\mathcal{L}}=j_{\text{surf}}\dot{\mathcal{M}},\end{array}$$(18)

where j_surf is the specific angular momentum of the distorted surface and $\dot{\mathcal{M}}$ the global mass-loss rate (Paxton et al. 2019). The local scheme computes angular momentum loss according to Eqs. (12) and (13), taking into account both the anisotropic wind and surface deformation caused by rotation.

We compute models with a chemical mixture representing the Large Magellanic Cloud as in Brott et al. (2011). Two initial masses of 10 M_⊙ and 20 M_⊙ are chosen to straddle the bi-stability jump. We run models with an initial equatorial rotation velocity of 300 km s⁻¹ until core hydrogen depletion, defined as a central hydrogen mass fraction of 1 × 10⁻⁴. The chosen rotation velocity represents the typical value for early B-stars found in the Large Magellanic Cloud (Dufton et al. 2013) and corresponds to an initial critical rotation fraction of around 0.45 for both masses.

It is also useful to assess numerical models with varying rotation rates at a fixed point in their evolution. To this end we run models with very small timesteps until the model has burnt 3% by mass of its initial supply of hydrogen in the core (X_c = 0.7169). This is approximately the earliest point at which the model finds itself in thermal and nuclear equilibrium and hence is a good point in the star’s evolution to investigate. We shall term the point when X_c = 0.7169 the zero-age-main-sequence.

Our models are numerically stable until initial critical velocity fractions, ω, of around 0.65, so to investigate stars with faster rotation, an extrapolation is performed. The local wind scheme requires as inputs the stellar mass, rotation rate, polar radius and luminosity.

The polar radius is assumed to be invariant to rotation, so this is known from a non-rotating model. For the luminosity, we extrapolate from the slower rotating models as described in Appendix B up to ω = 0.9. Using the four named quantities, the effective gravity and effective temperature profiles can be calculated (as outlined in Sect. 2.1.1), and resultingly the surface mass-flux. Thus we may investigate the behaviour of our scheme for very fast rotating stars on the zero-age-main-sequence despite not having stellar models at these rotations. Stars born with moderate rotation may evolve to rotate at high critical velocity fractions owing to internal structural changes (Hastings et al. 2020b), so it is important to check the behaviour of our scheme at near critical rotations.

3. Results

3.1. Mass loss on the zero-age-main-sequence

Figure 1 shows the surface mass flux as a function of colatitude for 10 M_⊙ and 20 M_⊙ models rotating at various rates. All models displayed have a central hydrogen mass fraction of 0.7169, equating to 97% of the initial hydrogen mass fraction. It is seen that for slow rotation, mass flux is stronger at the poles and weaker at the equator. This occurs because rotation results in a hotter pole, relative to the non-rotating case, and a cooler equator, and stellar winds are very sensitive to effective temperature changes.

Fig. 1. Mass-loss rate per unit surface area as a function of colatitude, θ, for 10 M⊙ (left panel) and 20 M⊙ models (right panel) at various rotation rates. All models have burnt 3% by mass of their initial hydrogen (i.e., Xc = 0.7169). Critical velocity fraction, ω, is depicted in the legend.

For faster rotating 10 M_⊙ models, mass-flux experiences a jump at colatitudes between 60 and 80° caused by the bi-stability jump. Moving from pole to equator across the stellar surface, effective gravity and hence effective temperature decrease. At some point, the effective temperature subceeds the ‘jump temperature’ at which Fe IV recombines to Fe III causing a sudden, dramatic increase in the mass flux, as evidenced in the left panel of Fig. 1. The 20 M_⊙ model does not undergo the same phenomenon as here the effective temperature always exceeds the jump temperature.

The global mass-loss rate depends on both the surface mass-flux and the stellar surface area. For the wind recipe of Vink et al. (2001) used in this work, provided the ionisation equilibrium does not change significantly, faster rotation is seen to cause a decreasing surface mass-flux at the equator. Also rotation increases the surface area of the equatorial region due to the equatorial bulge. These two effects can offset one another, causing the global mass-loss rate to be roughly independent of rotation, as exemplified by models shown in Fig. 2. We note that because of the relatively small area covered by the polar region, the polar surface mass-flux does not contribute significantly to the global mass-loss rate.

Fig. 2. Upper panels: global mass-loss rates as a function of critical velocity fraction, ω, for 10 M⊙ (left panel) and 20 M⊙ models (right panel). Predictions of the local scheme, where surface mass flux is determined by Eq. (16) and the global mass-loss rate given by Eq. (11), are given in red. The standard scheme, where the mass loss rates are increased by Eq. (17) is depicted in black. All models have burnt 3% by mass of their initial hydrogen (i.e., Xc = 0.7169). Dotted lines represent mass-loss rates calculated from extrapolation (see Sect. 2.3). Lower panels: ratios of rotating star mass-loss to non-rotating star mass loss, $\dot{\mathcal{M}}(\omega )/\dot{\mathcal{M}}(\omega =0)$, for each combination of mass-loss model and initial mass as given in each panel.

For the 10 M_⊙ model in the local scheme, mass loss decreases slightly with faster rotation, until the bi-stability jump comes into effect at ω ≈ 0.5 and drives mass-loss rate up. In contrast, the 20 M_⊙ model displays almost no change in mass-loss rates until ω ≈ 0.3 and then a small increase thereafter, due to the effect of the growing surface area of the equator dominating over the diminished equatorial surface mass-flux. We note that, except for models affected by the bi-stability jump, our local scheme produces slightly weaker winds than the standard scheme. Our models show that excluding the effects of the bistability jump, mass-loss rates of a rotating star are not predicted to be significantly different to that of an equivalent non-rotating star.

At very high initial rotation rates, our estimates of the mass-loss rate from extrapolation of the luminosity show that for the 20 M_⊙ model, the increase in mass-loss rate is modest, 10% at most. Whereas the cooler 10 M_⊙ model displays mass-loss rate enhancement of a factor 9 at ω = 0.9. We are thus confident that our scheme behaves reasonably at near-critical velocities.

3.2. Evolutionary models

Our evolutionary models are presented in Fig. 3 where panels a and b show the evolution of global mass-loss rates. In the standard scheme the bi-stability jump is implemented as a sharp jump, however in the local scheme the jump is more gradual, owing to the fact that as the star cools, the region of the surface affected by the jump grows, causing the mass-loss to also gradually increase.

Fig. 3. Upper panels: evolution of global mass-loss rate as a function of time. The blue dashed line shows the mass-loss rate computed from a non-rotating model. Central panels: evolution of normalised specific angular momentum loss, as given by Eq. (19) (see text for details). In the limit of slow rotation this quantity is equal to unity, shown by the blue dashed line. Lower panels: evolution of the equatorial rotational velocity as a function of time. The left panels show a 10 M⊙ model, the right panels a 20 M⊙ model. Predictions of the local scheme, where surface mass flux is determined by Eq. (16) and the global mass-loss rate given by Eq. (11), are given in red. The standard scheme, where the mass loss rates are increased by Eq. (17) is depicted in black. All rotating models have initial equatorial rotational velocities of 300 km s−1.

Panels c and d of Fig. 3 show the normalised specific angular momentum loss of our evolutionary models, given as

$$\begin{array}{c}\hfill \frac{\dot{\mathcal{L}}}{\frac{2}{3}\dot{\mathcal{M}}\mathrm{\Omega }{R}_{\text{eq}}^2}.\end{array}$$(19)

This is a unitless quantity that describes the strength of angular momentum loss independently of the rotation rate and mass-loss rate. A spherical star with an isotropic wind (i.e., a slowly rotating star) has a normalised specific angular momentum loss of 1. Values larger than unity imply that the star is losing more angular momentum per unit mass than the spherically symmetric case and that spin-down will occur more rapidly. This quantity is sensitive to both the anisotropy of the wind and the deformation of the star. We see that away from the bi-stability jump temperatures, models using the local scheme suffer lower angular momentum losses than the standard scheme. This reduced normalised specific angular momentum loss means that stars may be able to maintain faster rotation rates. The opposite is true when mass-flux across the stellar surface is increased due to the bi-stability jump, because the model loses large quantities of mass from the equatorial regions.

Panels e and f of Fig. 3 show the equatorial velocities of our models. Comparing the velocities near the end of the main-sequence, we see that the local scheme displays larger rotational velocities, due to the generally lower mass-loss and normalised specific angular momentum losses as shown in the upper two panels. The effect is greatest in the 20 M_⊙ model, with velocities increased by roughly 10% compared to the standard scheme.

4. Uncertainties

When attempting to describe a two-dimensional phenomenon with a one-dimensional model there are inevitably shortcomings. Most stellar evolution codes compute the structure of a rotating star by applying certain corrections to the stellar structure equations that are designed to produce the average properties along an isobar (for a detailed description see Heger et al. 2000). This approximation may break down under certain conditions, for example when the surface temperature at the equator is cool enough for helium-I to form yet the pole it is not, the opacity will vary greatly over the stellar surface causing different physical conditions at the equator and pole. In such a case, average quantities will not capture this diversity and may lead to different structures as computed by one and two-dimensional models.

A weakness of our wind scheme is that to determine the local mass-flux, we use a mass-loss recipe that was calculated for non-rotating stars. Such a recipe naturally ignores rotational phenomena like non-spherical geometry and the effects of limb-darkening. What is more, the ionisation of the wind is expected to be sensitive to radiation from various latitudes on the stellar surface (Petrenz & Puls 2000), which could effect the mass-loss rates.

A fundamental assumption of our scheme is that the wind is launched from the stellar surface and moves parallel to the photon flux (which is assumed to correspond to the direction of the effective gravity). In reality, a wind is continually accelerated until it reaches the terminal wind velocity and during this acceleration a wind particle may be influenced by photons streaming at an oblique angle to the stellar surface. This would introduce a non-radial line force (particularly in combination with a polar-angle dependent velocity field), which may alter the wind structure and angular momentum content (Owocki et al. 1996; Gayley & Owocki 2000). What is more, in our model the effective gravity, and hence flux, have a non-radial direction, while most one-dimensional mass-loss recipes assume a purely radial flux.

When running models at large critical velocity fractions, caution must be exercised, as one may be extrapolating from the non-rotating wind recipe. For any given wind recipe, there are bounds in which the input parameters are valid and it is entirely possible that under extreme rotation the local surface conditions fall outside of the prescribed bounds. Should this occur, a second suitable wind recipe could be used to give the local mass flux for the effected regions, for example a cool star wind recipe may be appropriate for describing the equatorial wind.

5. Applicability of the local mass-loss scheme

The pre-requisites for the mass-loss scheme presented here are that the star’s shape needs to be well described by the Roche potential and the gravity darkening law of Espinosa Lara & Rieutord (2011) must be valid. This is in general true for both convective and radiative stars that do not have near-Eddington luminosities and are rotating sub-critically (Espinosa Lara & Rieutord 2011), however there are further cases where these conditions are not met and other special cases which will be discussed here.

For very luminous stars, the radiative acceleration may facilitate the unbinding of material from the stellar surface at lower rotation velocities than the Keplerian velocity. As this effect is ignored in our formalism, our scheme is not appropriate for very luminous objects. In light of the findings of Maeder & Meynet (2000), we would conservatively advise the limit of applicability to be 60% of the Eddington luminosity. Improved gravity darkening laws and a more detailed account of the stellar surface opacity could possibly allow for applying our scheme also at higher Eddington factors, but this still needs to be investigated (see Sect. 2.2).

Furthermore, luminous stars may suffer the effects of inflation whereby radiation pressure “inflates” the star, producing a very tenuous, extended envelope (Ishii et al. 1999; Sanyal et al. 2015). If the radiation pressure deviates from spherical symmetry, then the strength of inflation will vary according to latitude, suggesting that the star’s shape is not well described by the Roche potential. Our mass-loss scheme is therefore not applicable to inflated stellar models. Inflation is expected to occur at masses above 30 M_⊙ for stars with galactic metallicity, but for much higher masses at lower metallicities (Sanyal et al. 2017).

A star may suffer the effects of additional forces which can alter the surface effective gravity beyond the Roche potential. Examples include radial pulsations, accelerations from rapid expansion or contraction and a close binary companion. Such cases would need to be dealt with separately, although our scheme could be extended to them.

While our scheme may be applied to rapid rotators, once a star reaches critical velocity, evolutionary models demand that the star lose enough angular momentum to maintain sub-critical rotation. It is not entirely clear how this may happen, there are several possibilities. The star may undergo a “mechanical mass-loss episode”, losing the required angular momentum through increased mass loss at the equator only (Granada et al. 2013). The other extreme is to lose angular momentum via an isotropic wind, as is currently done in MESA models, but one may also prescribe for mass to be lost from the surface in any configuration. For at least some fast rotating stars in nature, a circumstellar decretion disc forms that can efficiently drain angular momentum from the star (Krtička et al. 2011). For lower-mass stars (M ⪅ 10 M_⊙), the required angular momentum loss rates, and correspondingly required mass loss rates are low (Granada et al. 2013), hence the evolution of the star is largely insensitive to the mechanics of angular momentum loss at the critical velocity. This is not true for more massive stars, which can lose upwards of 10% of their initial mass from rotating critically (cf. Table 1 of Granada et al. 2013), so how exactly angular momentum is drained from a critical rotator becomes important. Therefore we advise caution when stellar models achieve critical rotation.

A crucial aspect of our formulation is that it demands that the wind is sensitive only to the local conditions of where on the surface it was launched from. An example where this condition is violated is the dust-driven winds of asymptotic giant branch stars. Global pulsations may lead to dust formation in the outer atmosphere, which is essential for the wind driving (Winters et al. 2000).

6. Comparison to other studies

Several authors have investigated the problem of stellar winds and rotation by directly taking into account rotation specific physics. The prescription presented here is better described as an adaptation of a wind model for non-rotating stars, so it is useful to compare our results to previous studies.

It has been reported that radiation-driven winds are most strongly affected by gravity darkening directly beneath the point from which the wind was launched (Cranmer & Owocki 1995). This suggests that the wind is only sensitive to the point form which it is launched, justifying our use of a non-rotating wind model as our basis. It is also encouraging as limb-darkening, which is not accounted for in our prescription, is deemed unimportant (Cranmer & Owocki 1995).

Petrenz & Puls (2000) calculated wind models using the concept of a mean irradiating atmosphere and found the winds to have a prolate structure, with increased mass-flux at the pole. Furthermore, for B-type stars, rotation is predicted to diminish mass-loss rates, with models rotating at around 80% of critical velocity displaying mass-loss rates a few percent lower than corresponding non-rotating models (cf. Table 4 of Petrenz & Puls 2000). Similarly, we predict a very weak rotation dependence on mass-loss rates (away from the bi-stability jump), although our models can show enhanced or reduced winds depending on the stellar parameters. Müller & Vink (2014) also find that mass-loss actually diminishes due to the effects of rotation, in contrast to the rotationally enhanced wind schemes.

Pelupessy et al. (2000) focused on B[e] stars using models including the bi-stability jump effect. They report that rotation, in general, enhances mass-flux from the poles but hardly changes that of the equator. The spatial variation in mass-flux predicted by Pelupessy et al. (2000) shows a discontinuity owing to the bi-stability jump, albeit not as steep as in our results (cf. Fig. 9 with our Fig. 1). Pelupessy et al. (2000) predict the winds of a 20 M_⊙ star to grow stronger with rotation, with rotation at 60% of critical velocity boosting mass-loss by 16% compared to the non-rotating case (cf. Table 3).

The works mentioned above computed only stationary models, however different stellar parameters were used in each case. For example the 20 M_⊙ model of Petrenz & Puls (2000) had a radius of 20 R_⊙, while that of Pelupessy et al. (2000) was more than twice as large, 47 R_⊙. The fact that the resulting relationships between rotation and mass-loss rates disagree is therefore not surprising.

The work of Gagnier et al. (2019a) differs from this study twofold. Firstly, the two-dimensional ESTER code (Espinosa Lara & Rieutord 2013) was used to compute the stellar structure, whereas here we rely on a one-dimensional code. Secondly, the local mass-flux was calculated by calibrating the one-dimensional CAK theory (Castor et al. 1975) to the wind recipe of Vink et al. (2001). Given that models presented here and in Gagnier et al. (2019a) are based on the same wind recipe of Vink et al. (2001), a comparison between the two will highlight differences in the underlying methods.

To compare our scheme to that of Gagnier et al. (2019a), a 15 M_⊙ model has been computed with solar metallicity as in Brott et al. (2011) and initial critical angular velocity fraction, ω = 0.5. Figure 4 compares the results of the two methods. Firstly, we see that the equatorial effective temperatures predicted by both models agree to within around 2000 K, and to within several hundred Kelvin during the early evolution. This discrepancy can largely be credited to the differences of the structures predicted by two and one-dimensional models and the implementation of rotational mixing.

Fig. 4. Comparison of the global mass-loss rates (left panel) and equatorial effective temperatures (right panel) predicted by this work and that of Gagnier et al. (2019a). Shown are the results of 15 M⊙ models with solar metallicity and initial critical velocity fraction, ω = 0.5. Predictions of this work are plotted as grey lines, and those of Gagnier et al. (2019a) as cyan lines. The x-axis depicts the central hydrogen mass fraction normalised to the initial value. The panels in this figure are directly comparable to Figs. 13 and 14 of Gagnier et al. (2019a).

A disadvantage of the method of Gagnier et al. (2019a) is that the required calibration of the surface mass-flux is sensitive to the strength of surface gravity, meaning that properly, (as stars evolve to lower surface gravities) a new calibration must be made at every timestep (see Sect. 4.2 of Gagnier et al. 2019b). However, as the calibration is onerous, it was only carried out for models on the zero-age-main-sequence, meaning that “the local mass-flux may be underestimated by a factor ∼1.7 at most”. Our scheme does not suffer from this issue, which may explain partly why our model predicts slightly higher mass-loss rates in panel a of Fig. 4. Both models however show the same general trend, with a gradual increase in mass-loss rates once the equatorial effective temperature cools below 22–23 kK. The relative increase in mass-loss rates brought about by the bistability jump is approximately the same in both models. The jump temperature differs slightly in the two models because the jump temperature is sensitive to the stellar luminosity (see Eq. (16) of Gagnier et al. 2019a). Both models naturally have different luminosities owing to their different structures, as mentioned earlier. We find it encouraging that our relatively simple scheme behaves similarly to a two-dimensional, more advanced model.

7. Conclusions

We have presented a new and simple to implement prescription for the mass and angular momentum loss rates of rotating massive stars. This represents an improvement over the widely used rotationally enhanced mass-loss schemes as we calculate the two-dimensional mass-flux over the stellar surface and are able to compute the angular momentum loss resulting from an anisotropic wind originating from a distorted star. Our method involves using a mass-loss recipe for non-rotating stars to determine the local mass-flux across the surface of a rotating star, which is then integrated to give global mass and angular momentum loss rates.

In general we notice that, away from the bi-stability jump temperature, mass-loss rates are slightly diminished compared to the rotationally enhanced mass-loss scheme. The local mass-flux scheme has the effect of smoothing out the bi-stability jump as the increase in mass-loss rate is implemented locally on the star’s surface, not globally, also observed in the models of Gagnier et al. (2019a). Our models show that the presence of the bistability-jump causes a strong relationship between rotation and mass-loss. If the bistability-jump does not in fact operate in nature, as suggested by theoretical wind models of Björklund et al. (2023), moderate and even fast rotation is not predicted to strongly alter mass-loss compared to the non-rotating case. We see evidence that the detailed relationship between rotation and wind strength is complex, with mass-loss rates being either decreased or increased depending on the surface properties of the star.

Our aim to provide a scheme for one-dimensional stellar evolution codes will of course mean that simplifications must be made. In spite of our scheme’s shortcomings, comparisons with similar, more physically comprehensive works deliver a broad agreement in global mass-loss rates.

Our methods are relevant to several areas of stellar astrophysics where the evolution of angular momentum plays a decisive role. For example, Be stars are known to be fast rotators, so are expected to have strongly anisotropic winds and large distortions. Models such as those presented here may therefore be used to investigate evolutionary properties of Be stars. Secondly, the evolution of models along a chemically homogeneous pathway can be interrupted by spin-down caused by stellar winds (Yoon et al. 2006). The calculations presented here suggest that angular momentum loss has been generally overestimated in stellar models, suggesting that chemically homogeneous evolution (e.g., Hastings et al. 2020a) may be more common or easier to maintain than previously thought. In the future, our wind prescription may be implemented in the next grids of stellar evolution models in order to gain further insights into the physics of rotating massive stars.

Acknowledgments

The authors are grateful to Pablo Marchant for providing the MESA input files that form the basis of this work. Thanks also go to the anonymous referee who helped to improve this manuscript.

References

Appendix A: Shapes of rotating stars

Here we derive the shape of the surface of a rotating star. We assume that the star is well described as a point mass enclosed by a massless envelope and that the polar radius in unaffected by rotation. The surface of the star is then an equipotential given by the Roche potential and reads

$$\begin{array}{c}\hfill \frac{\mathit{GM}}{r(\theta )}+\frac{1}{2}{\mathrm{\Omega }}^{2}r{(\theta )}^2{\text{sin}}^2(\theta )=\text{constant},\end{array}$$(A.1)

with M representing the stellar mass, Ω the angular velocity (which is assumed to be constant across both the surface and through the interior of the star), r(θ) the radial co-ordinate and θ the co-latitude. We parametrise the strength of rotation with the Keplerian angular velocity, defined using the equatorial radius R_e as

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\text{Kep}}^2=\frac{\mathit{GM}}{R_e^3}.\end{array}$$(A.2)

We let

$$\begin{array}{c}\hfill \omega =\mathrm{\Omega }/{\mathrm{\Omega }}_{\text{Kep}}.\end{array}$$(A.3)

It is important to stress that the derivation that follows is only valid for the above parametrisation of rotation.

Combining Eqns. A.1, A.2 and A.3 gives

$$\begin{array}{c}\hfill 1+\frac{1}{2}{\omega }^2\frac{r^3}{R_e^3}{\text{sin}}^2(\theta )=\frac{r}{R_e}(1+\frac{{\omega }^2}{2}),\end{array}$$(A.4)

which after defining

$$\begin{array}{c}\hfill \stackrel{\sim }{r}=r/R_e,\end{array}$$(A.5)

further simplifies to

$$\begin{array}{cc}\hfill {\stackrel{\sim }{r}}^3-\frac{2+{\omega }^2}{{\omega }^2{\sin }^2(\theta )}\stackrel{\sim }{r}+\frac{2}{{\omega }^2{\sin }^2(\theta )}& =0.\hfill \end{array}$$(A.6)

Eq. A.6 is a cubic in the form x³ + px + q = 0 (known as a depressed cubic) and has the general solution

$$\begin{array}{c}\hfill x_k=2\sqrt{\frac{-p}{3}}\text{cos}[\frac{1}{3}\text{arccos}\left(\frac{3q}{2p}\sqrt{\frac{-3}{p}}\right)-\frac{2\pi k}{3}]\end{array}$$(A.7)

for k = 0, 1, 2 corresponding to the 3 cubic roots (Zwillinger 1996).

Here we have $p=-\frac{2+{\omega }^2}{{\omega }^2{\sin }^2(\theta )}$ and $q=\frac{2}{{\omega }^{2}si{n}^2(\theta )}$ giving

$$\begin{array}{c}\hfill {\stackrel{\sim }{r}}_k=2\sqrt{\frac{2+{\omega }^2}{3{\omega }^2{\text{sin}}^2(\theta )}}\text{cos}[\frac{1}{3}\text{arccos}(-\frac{3}{2+{\omega }^2}\sqrt{\frac{3{\omega }^2{\text{sin}}^2(\theta )}{2+{\omega }^2}})-\frac{2\pi k}{3}].\end{array}$$(A.8)

The periodicity of the cosine function means that if y = cos(x) then −y = cos(x + nπ) where n is an integer so that arccos(−y) = x + nπ and x = arccos(y)−nπ. Thus Eq. A.8 becomes

$$\begin{array}{c}\hfill {\stackrel{\sim }{r}}_k=2\sqrt{\frac{2+{\omega }^2}{3{\omega }^2{\sin }^2(\theta )}}\text{cos}[\frac{1}{3}\text{arccos}\left(\frac{3}{2+{\omega }^2}\sqrt{\frac{3{\omega }^2{\sin }^2(\theta )}{2+{\omega }^2}}\right)-\frac{2\pi k}{3}-n\pi ].\end{array}$$(A.9)

To be physical, the solution must be independent of n, which is only achieved when k = 1. Choosing n = −1 results in

$$\begin{array}{c}\hfill {\stackrel{\sim }{r}}_{k=1}=2\sqrt{\frac{2+{\omega }^2}{3{\omega }^2{\sin }^2(\theta )}}\text{cos}[\frac{1}{3}\text{arccos}\left(\frac{3}{2+{\omega }^2}\sqrt{\frac{3{\omega }^2{\sin }^2(\theta )}{2+{\omega }^2}}\right)+\frac{\pi }{3}].\end{array}$$(A.10)

From Eqn. A.4 one can deduce the ratio of equatorial to polar radii, $\frac{R_{\mathrm{e}}}{R_p}$ as

$$\begin{array}{c}\hfill \frac{R_e}{R_p}=1+\frac{1}{2}{\omega }^2.\end{array}$$(A.11)

The final solution is arrived at by combining Eqns. A.10 and A.11 and reads

$$\begin{array}{c}\hfill \frac{r(\omega ,\theta )}{R_{\text{p}}}=(2+{\omega }^2)\sqrt{\frac{2+{\omega }^2}{3{\omega }^2{\sin }^2(\theta )}}\text{cos}[\frac{1}{3}\text{arccos}\left(\frac{3}{2+{\omega }^2}\sqrt{\frac{3{\omega }^2{\sin }^2(\theta )}{2+{\omega }^2}}\right)+\frac{\pi }{3}].\end{array}$$(A.12)

As expected, the expression above gives the equatorial radius at critical rotation to be 1.5 times the polar radius ($\frac{r(\omega =1,\theta =\pi /2)}{R_p}=1.5$). We note that Eq. A.12 differs from Equation 26 used by Cranmer & Owocki (1995) owing to the use of different definitions of the critical velocity Cranmer & Owocki (1995) use ${\mathrm{\Omega }}_{\text{crit}}^2\propto M/{(1.5R_p)}^3$.

Appendix B: Surface properties of fast rotators close to the zero-age-main-sequence

Our numerical models are unable to compute the structure of a rotating star with initial rotation exceeding around 70% of the critical velocity. However the wind properties of very fast rotating stars on the zero-age-main-sequence may be determined via extrapolation. The wind scheme presented in this paper requires knowledge of the effective gravity and effective temperature profile of a star. To calculate these profiles, only the polar radius, luminosity, mass and rotation rate are required.

To calculate the luminosity of fast rotators on the zero-age-main-sequence, we extrapolate linearly from our models with initial critical fractions between 0.4 and 0.6. We extrapolate the luminosity normalised to the value of the non-rotating model against the initial fraction of critical velocity, where all values are defined at the point when the central hydrogen mass fraction decreases by 3% from its initial value (this we term the zero-age-main-sequence, as it is the earliest point in which the models find themselves in equilibrium). This extrapolation is shown in Fig. B.1 for 10 and 20M_⊙ models. We find that rotation rate and luminosity are inversely proportional. This effect is rather weak, with luminosity decreasing by approximately 6% at 60% of critical rotation. The extrapolations suggest that at up to 90% of the critical velocity, the luminosity is reduced by no more than 10% compared to the non-rotating case.

Fig. B.1. Variation of luminosity (black line), L, and polar radius (yellow line), Rp, normalised to the values of a non-rotating star as a function of critical velocity fraction, ω. The left panel shows models with masses 10M⊙ and the right panel 20M⊙. All models have burnt 3% by mass of their initial hydrogen (i.e. Xc = 0.7169). Each cross represents a value computed by a MESA model. Dashed lines show extrapolated values.

In line with the assumptions of the Roche potential, the polar radius is not affected by rotation, so this value may be assumed from a non-rotating MESA model, depicted graphically by the horizontal dashed yellow line in Fig. B.1. Indeed, as evidenced by Fig. B.1 the numerical models predict that to within a few percent, the polar radius remains unchanged by rotation.

The results of our extrapolations compare favourably to the computed stellar structures of Ekström et al. (2008). One-dimensional models of a 20M_⊙ star predict that the luminosity decreases by around 8% over the course of being spun up from stationary to critical rotation (Ekström et al. 2008, Fig. 5). While the polar radius is judged to shrink very slightly with increasing rotation, but this is at most a 2% effect (Ekström et al. 2008, Fig. 2), thus justifying the assumption of the Roche potential.

All Figures

	Fig. 1. Mass-loss rate per unit surface area as a function of colatitude, θ, for 10 M⊙ (left panel) and 20 M⊙ models (right panel) at various rotation rates. All models have burnt 3% by mass of their initial hydrogen (i.e., Xc = 0.7169). Critical velocity fraction, ω, is depicted in the legend.
In the text

Fig. 2. Upper panels: global mass-loss rates as a function of critical velocity fraction, ω, for 10 M⊙ (left panel) and 20 M⊙ models (right panel). Predictions of the local scheme, where surface mass flux is determined by Eq. (16) and the global mass-loss rate given by Eq. (11), are given in red. The standard scheme, where the mass loss rates are increased by Eq. (17) is depicted in black. All models have burnt 3% by mass of their initial hydrogen (i.e., Xc = 0.7169). Dotted lines represent mass-loss rates calculated from extrapolation (see Sect. 2.3). Lower panels: ratios of rotating star mass-loss to non-rotating star mass loss, $\dot{\mathcal{M}}(\omega )/\dot{\mathcal{M}}(\omega =0)$, for each combination of mass-loss model and initial mass as given in each panel.

In the text

Fig. 3. Upper panels: evolution of global mass-loss rate as a function of time. The blue dashed line shows the mass-loss rate computed from a non-rotating model. Central panels: evolution of normalised specific angular momentum loss, as given by Eq. (19) (see text for details). In the limit of slow rotation this quantity is equal to unity, shown by the blue dashed line. Lower panels: evolution of the equatorial rotational velocity as a function of time. The left panels show a 10 M⊙ model, the right panels a 20 M⊙ model. Predictions of the local scheme, where surface mass flux is determined by Eq. (16) and the global mass-loss rate given by Eq. (11), are given in red. The standard scheme, where the mass loss rates are increased by Eq. (17) is depicted in black. All rotating models have initial equatorial rotational velocities of 300 km s−1.

In the text

Fig. 4. Comparison of the global mass-loss rates (left panel) and equatorial effective temperatures (right panel) predicted by this work and that of Gagnier et al. (2019a). Shown are the results of 15 M⊙ models with solar metallicity and initial critical velocity fraction, ω = 0.5. Predictions of this work are plotted as grey lines, and those of Gagnier et al. (2019a) as cyan lines. The x-axis depicts the central hydrogen mass fraction normalised to the initial value. The panels in this figure are directly comparable to Figs. 13 and 14 of Gagnier et al. (2019a).

In the text

Fig. B.1. Variation of luminosity (black line), L, and polar radius (yellow line), Rp, normalised to the values of a non-rotating star as a function of critical velocity fraction, ω. The left panel shows models with masses 10M⊙ and the right panel 20M⊙. All models have burnt 3% by mass of their initial hydrogen (i.e. Xc = 0.7169). Each cross represents a value computed by a MESA model. Dashed lines show extrapolated values.

In the text