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

TL;DR

This paper investigates the critical angular velocity and mass loss dependence on rotation in massive stars with radiation-driven winds, using 2D models to understand their evolution and wind regimes.

Contribution

It introduces a combined analysis of radiative acceleration effects and a local surface mass flux model for rotating stars, highlighting the impact on mass loss regimes and critical velocity.

Findings

Radiative acceleration mildly affects critical angular velocity for stars under 40 solar masses.

Mass loss is significantly enhanced in two-wind regimes around the equator.

Detecting the wind regime transition observationally is challenging due to viewing angle effects.

Abstract

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.

Tables

Table 1. Table 1: Summary of the main results for ESTER 2D models of a 15 ​ M ⊙ 15 subscript 𝑀 direct-product 15\leavevmode\nobreak\ M_{\odot} star with Z = 0.02 𝑍 0.02 Z=0.02 at ZAMS. The first column reports the ratio of equatorial angular velocity to Keplerian angular velocity ω 𝜔 \omega , the second column reports the critical angular velocity ratio Ω e ​ q / Ω c subscript Ω 𝑒 𝑞 subscript Ω 𝑐 \Omega_{eq}/\Omega_{c} , the third column is the global mass-loss rate log ⁡ M , ˙ \log\dot{M,} and the last column gives the ratio between total angular momentum and angular momentum loss rate log ⁡ T ℒ = log ⁡ ℒ / ℒ ˙ subscript T ℒ ℒ ˙ ℒ \log{\rm T}_{\mathcal{L}}=\log\mathcal{L}/\dot{\mathcal{L}} . M ˙ ˙ 𝑀 \dot{M} is in M ⊙ ⋅ yr − 1 ⋅ subscript 𝑀 direct-product superscript yr 1 M_{\odot}\cdot{\rm yr}^{-1} and T ℒ subscript T ℒ {\rm T}_{\mathcal{L}} in yr.

$\omega$	${\Omega }_{eq}/{\Omega }_c$	$\log \dot{M}$	$\log T_L$
0	0	$-8.40$	–
0.1	0.105	$-8.41$	8.60
0.2	$0.210$	$-8.42$	$8.61$
0.3	$0.313$	$-8.43$	$8.62$
0.4	$0.415$	$-8.44$	$8.63$
0.5	$0.516$	$-8.46$	$8.64$
0.6	$0.615$	$-8.47$	$8.65$
0.7	$0.714$	$-8.48$	$8.65$
0.8	$0.812$	$-8.49$	$8.64$
0.9a	$0.909$	$-8.27$	$8.35$
0.95a	$0.957$	$-8.25$	$8.22$
0.99a	$0.995$	$-8.19$	$8.13$

Equations

$${\mathbf{g}}_{\mathrm{tot}}={\mathbf{g}}_{\mathrm{eff}}+{\mathbf{g}}_{\mathrm{rad}},$$

$${\mathbf{g}}_{\mathrm{rad}}=\frac{\kappa \mathbf{F}}{c},$$

$${\Omega }_k=\sqrt[]{\frac{GM}{R_{eq}^3}},$$

$$\Gamma =\frac{\kappa L}{4\pi cGM}$$

$$\mathbf{F}=-\rho \chi \frac{dT}{dP}{\mathbf{g}}_{\mathrm{eff}},$$

$$\rho \chi \frac{dT}{dP}=\frac{L(P)}{4\pi G{M}_{\star }},$$

$$M_{\star }=M\left(1-\frac{{\Omega }^2}{2\pi G{\rho }_{\mathrm{m}}}\right),$$

$$\mathbf{F}=-\frac{L(P)}{4\pi G{M}_{\star }}{\mathbf{g}}_{\mathrm{eff}}.$$

$$\mathbf{F}=-\frac{L(P)}{4\pi G{M}_{\star }}(1-\zeta (\theta )){\mathbf{g}}_{\mathrm{eff}}.$$

$$F(\theta )\sim g_{\mathrm{eff}}(\theta {)}^{4\beta },$$

$$\mathbf{F}=-f(r,\theta ){\mathbf{g}}_{\mathrm{eff}},$$

$$\underset{r\to 0}{\lim }f(r,\theta )=\frac{L}{4\pi GM}.$$

$$f(r,\theta )=\frac{L}{4\pi GM}\frac{{\tan }^2\psi (r,\theta )}{{\tan }^2\theta },$$

$$\cos \psi +\ln \tan (\psi /2)=\frac{1}{3}{\omega }^2{r}^3{\cos }^3\theta +\cos \theta +\ln \tan (\theta /2).$$

$$\omega =\frac{\Omega }{{\Omega }_k}$$

$$f(r=1,\pi /2)=\frac{L}{4\pi GM}{\left(1-\frac{{\Omega }^2{R}_{eq}^3}{GM}\right)}^{-2/3},$$

$$\mathbf{F}(R_{eq},\pi /2)=-\frac{L}{4\pi GM}{\left(1-{\omega }^2\right)}^{-2/3}{\mathbf{g}}_{\mathrm{eff}}.$$

$$\mathbf{F}(R_{eq},\pi /2)\simeq -\frac{L}{4\pi GM\left(1-\frac{2}{3}{\omega }^2\right)}{\mathbf{g}}_{\mathrm{eff}}.$$

$${\mathbf{g}}_{\mathrm{tot}}={\mathbf{g}}_{\mathrm{eff}}+{\mathbf{g}}_{\mathrm{rad}}=0$$

$${\mathbf{F}}_{\lim }=-\frac{c}{\kappa }{\mathbf{g}}_{\mathrm{eff}}.$$

$${\Gamma }_{\Omega }(\theta )=\frac{F(\theta )}{F_{\lim }(\theta )}=\frac{\kappa (\theta )}{c}f(r=1,\theta ).$$

$${\mathbf{g}}_{\mathrm{tot}}={\mathbf{g}}_{\mathrm{eff}}\left[1-{\Gamma }_{\Omega }(\theta )\right].$$

$$g_{\mathrm{tot}}(\pi /2)=g_{\mathrm{eff}}(\pi /2)+g_{\mathrm{rad}}(\pi /2),$$

$$g_{\mathrm{rad}}(\pi /2)=-\frac{\kappa (\pi /2)L}{4\pi cGM}{\left(1-{\omega }^2\right)}^{-2/3}{g}_{\mathrm{eff}}(\pi /2),$$

$$g_{\mathrm{eff}}(\pi /2)=-R_{eq}{\Omega }_k^2\left(1-{\omega }^2\right).$$

$${\tilde{g}}_{\mathrm{tot}}(\pi /2)={\omega }^2+{\Gamma }_{eq}(1-{\omega }^2{)}^{1/3}-1,$$

$${\Gamma }_{\Omega }(\pi /2)=\frac{\kappa (\pi /2)L}{4\pi cGM}{\left(1-\frac{{\Omega }_c^2}{{\Omega }_k^2}\right)}^{-2/3}=1,$$

$${\Omega }_c={\Omega }_k\sqrt{1-{\Gamma }_{eq}^{3/2}}.$$

$$\frac{\delta F}{F}=\frac{|F_{\mathrm{ESTER}}-F_{\omega }|}{F_{\mathrm{ESTER}}}$$

$$\dot{M}=4\pi r^2\rho v=\mathrm{Cst},$$

Full text

*Aims. *The aim of this paper is to investigate two questions that are prerequisite to the study of the evolution of massive rapidly rotating stars: (i) What is the critical angular velocity of a star when radiative acceleration is significant in its atmosphere? (ii) How do mass and angular momentum loss depend on the rotation rate?

*Methods. *To investigate fast rotation, which makes stars oblate, we used the 2D ESTER models and a simplified approach, the $\omega$-model, which gives the latitudinal dependence of the radiative flux in a centrifugally flattened radiative envelope.

*Results. *We find that radiative acceleration only mildly influences the critical angular velocity, at least for stars with masses lower than 40\leavevmode\nobreak\ \mbox{\rm M{}_{\odot}}. For instance, a 15 $M_{\odot}$ star on the zero-age main sequence (ZAMS) would reach criticality at a rotation rate equal to 0.997 the Keplerian equatorial rotation rate. We explain this mild reduction of the critical angular velocity compared to the classical Keplerian angular velocity by the combined effects of gravity darkening and a reduced equatorial opacity that is due to the centrifugal acceleration. To answer the second question, we first devised a model of the local surface mass flux, which we calibrated with previously developed 1D models. The discontinuity (the 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.

1 Introduction

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

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

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

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

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

2 Critical angular velocity and the $\Omega\Gamma$-limit

2.1 The $\Omega\Gamma$-limit question

2.1.1 Some context

In stars more massive than $\sim 7\leavevmode\nobreak\ M_{\odot}$ that are close to solar metallicity, radiative acceleration plays a significant role in the (assumed) hydrostatic equilibrium. Total gravity is usually introduced,

[TABLE]

where effective gravity (or gravito–centrifugal acceleration) ${\boldsymbol{g}}_{\rm eff}$ is supplemented by radiative acceleration,

[TABLE]

where $\kappa$ is the flux-weighted opacity111Strictly speaking, $\kappa$ is the mass absorption coefficient, $\kappa=\mu/\rho$, with opacity $\mu=\lambda^{-1}$, where $\lambda$ is the mean free path of photons, and density $\rho$. However, in the following, we use, as is frequently done in the current context, the term “opacity”., which we approximate with the total Rosseland mean opacity, ${\boldsymbol{F}}$ the radiative flux, and $c$ the speed of light.

The so-called $\Omega\Gamma$-limit introduced by Maeder & Meynet (2000) is reached when ${\boldsymbol{g}}_{\rm tot}=\boldsymbol{0}$ somewhere on the stellar surface. It is associated with an actual critical angular velocity $\Omega_{c}$ that is different from the Keplerian angular velocity

[TABLE]

which is the break-up limit (or the $\Omega$-limit) when radiative acceleration can be neglected in total gravity. In this equation, $G$ is the gravitation constant, $M$ is the stellar mass, and $R_{eq}$ is its equatorial radius.

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

[TABLE]

where $L$ is the stellar luminosity, has no effect on the critical rotation because of gravity darkening. In an attempt to clarify the debate, Maeder (1999) and Maeder & Meynet (2000) (hereafter referred to as MMM) re-investigated the question and found two roots to the equation ${\boldsymbol{g}}_{\rm tot}=\boldsymbol{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 $\Omega_{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 ${\boldsymbol{F}}$ at some colatitude $\theta$ on the surface of a rotating star is proportional to the local effective gravity ${\boldsymbol{g}}_{\rm eff}$. For barotropic stars, this leads to

[TABLE]

where $\chi=4acT^{3}/(3\kappa\rho)$ is the radiative conductivity. Additionally, assuming solid-body rotation, Zahn (1992) obtained

[TABLE]

where

[TABLE]

and $\rho_{\rm 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

[TABLE]

In the case of shellular rotation, $\Omega\simeq\Omega(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

[TABLE]

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 $\zeta(\theta)\sim 10^{-2}$ so that according to this model, the ratio $F/g_{\rm eff}$ depends only mildly on colatitude.

2.1.2 Interferometric observations

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

[TABLE]

then $\beta<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 $\beta\simeq 0.25$ at low rotation rates, but that $\beta$ decreases to 0.13 when rotation approaches criticality. This behaviour has implications on the $\Omega\Gamma$-limit introduced by Maeder & Meynet (2000). These limitations prompt us to revisit this limit with ESTER 2D models.

2.2 The $\omega$-model

Before using full ESTER 2D models, it is worth considering the problem in light of the simplified $\omega$-model of Espinosa Lara & Rieutord (2011). The general purpose of the $\omega$-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

[TABLE]

where $f(r,\theta)$ is some function of the position to be determined. In this assumption, ${\boldsymbol{F}}$ and ${\boldsymbol{g}}_{\rm 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,\theta)$ can be determined from flux conservation equation, namely $\boldsymbol{\nabla}\cdot{{\boldsymbol{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

[TABLE]

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

[TABLE]

where $\psi(r,\theta)$ is obtained by solving

[TABLE]

In this equation $r$ has been scaled with the equatorial radius $R_{eq}$ and

[TABLE]

where $\Omega$ 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,

[TABLE]

so that the equatorial radiative flux reads

[TABLE]

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

[TABLE]

In this limit, where $R_{eq}\approx 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 $\omega$-model) expression (17) shows that the ratio $F/g_{\rm eff}$ diverges when the $\Omega$-limit $\omega=1$ is approached, while in its slow rotation approximation (18), $F/g_{\rm 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 $\Omega\Gamma$-limit.

2.3 Critical angular velocity at the $\Omega\Gamma$-limit: Ideas from the $\omega$-model

2.3.1 Preliminaries

With the expression of the radiative flux from the $\omega$-model, we can derive the critical angular velocity $\Omega_{c}$ corresponding to the $\Omega\Gamma$-limit. When this limit is reached, then

[TABLE]

somewhere at the surface of the star. As in MMM, we introduce a limiting flux from Eq. (2) and the $\Omega\Gamma$-limit condition ${\boldsymbol{g}}_{\rm tot}=\boldsymbol{0}$, namely

[TABLE]

From this expression, we define the rotation-dependent Eddington parameter $\Gamma_{\Omega}(\theta)$ as the ratio of the actual flux $F(\theta)$ obtained with the $\omega$-model, and the limiting flux, namely

[TABLE]

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

[TABLE]

The critical angular velocity $\Omega_{c}$ is reached if somewhere on the stellar surface ${\boldsymbol{g}}_{\rm tot}=\boldsymbol{0}$, that is, if there is a colatitude where either $\Gamma_{\Omega}(\theta)=1$ or ${\boldsymbol{g}}_{\rm eff}(\theta)=\boldsymbol{0}$.

2.3.2 Critical latitude: the equator

In all 2D models both effective temperature and effective gravity are minimum at the equator (e.g. Espinosa Lara & Rieutord, 2013). The solution ${\boldsymbol{g}}_{\rm eff}(\theta)=\boldsymbol{0}$ is therefore always reached first at the equator. We now focus on the $\Gamma_{\Omega}(\theta)=1$ solution.

We first observe that $\Gamma_{\Omega}(\theta)\propto\kappa(\theta)f(r=1,\theta)$ should be an increasing function of co-latitude, at least for (very) rapidly rotating stars. $f(r=1,\theta)$ indeed always increases with $\theta$ 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 $\kappa$ increases or decreases with decreasing $T_{\rm eff}$ (see appendix A for an attempt). Still, we expect the opacity to vary on the stellar surface, but much less than $f(r,\theta)$ near Keplerian angular velocity. According to the $\omega$-model, and because of the equatorial singularity, it is therefore very likely that the solution $\Gamma_{\Omega}=1$ is always first reached at the equator.

Maeder & Meynet (2000) came to the same conclusion regarding the location of the $\Gamma_{\Omega}(\theta)=1$ solution on the surface. However, they traced it back to an opacity effect, assuming that the latter increases with decreasing $T_{\rm eff}$ and thus is highest at the equator. According to ESTER 2D models, this is not the case for rotating stars (see below).

2.3.3 Unique critical angular velocity

Equation (22) shows that there are two solutions for ${\boldsymbol{g}}_{\rm tot}=\boldsymbol{0}$, and thus two possible critical angular velocities. However, as we show now, the $\omega$-model removes the ${\boldsymbol{g}}_{\rm eff}=\boldsymbol{0}$ root for the $\Omega\Gamma$-limit and thus points to a single critical angular velocity. According to the $\omega$-model at equator, Eq. (1) reads

[TABLE]

where

[TABLE]

and

[TABLE]

Here $g_{\rm eff}$ and $g_{\rm rad}$ are the radial components of the accelerations (thus positive when outwards). We can then write the equatorial total gravity scaled with $\Omega_{k}^{2}R_{eq}$ as

[TABLE]

where $\Gamma_{eq}$ is the standard Eddington parameter evaluated at the equator. From Eq. (24), we see that at the equator, the ratio $g_{\rm rad}/g_{\rm eff}$ increases as $(1-\omega^{2})^{-2/3}$ with increasing $\omega$, which also implies that if $g_{\rm eff}$ approaches [math] when $\omega\rightarrow 1$, $g_{\rm rad}$ will also tend to [math] but more slowly. Figure 1 shows the scaled total gravity, effective gravity, and radiative acceleration at the equator as a function of $\omega$ with $\Gamma_{eq}=0.5$. The total gravity at the equator has two zeros; the first root corresponds to $\Gamma_{\Omega}(\pi/2)=1,$ and the second root gives $g_{\rm eff}(\pi/2)=0$.

For sub-critical rotation (i.e. $g_{\rm tot}<0$ or $|g_{\rm eff}|>g_{\rm rad}$), the star is gravitationally bound. When we increase $\omega$, the equatorial effective gravity $|g_{\rm eff}|$ decreases faster than $g_{\rm rad}$, to the point where $|g_{\rm eff}|=g_{\rm rad}$ (equivalently, $\Gamma_{\Omega}(\pi/2)=1$), at this point, $\Omega=\Omega_{c}$ and $g_{\rm tot}=0$. Increasing $\omega$ even more would result in a radiative acceleration that surpasses the effective gravity at equator. When this happens, $g_{\rm tot}>0$ and the star is no longer gravitationally bound up to $\omega=1$ where the second root is reached. The solution $\Gamma_{\Omega}(\pi/2)=1$ is therefore always reached before $\omega=1$, when evolution (say) drives the growth of $\omega$.

2.3.4 Critical rotation given by the $\omega$-model

The main difference between the MMM model and ours, in addition to our unique critical angular velocity, comes from the latitudinal variation of $\Gamma_{\Omega}$. In MMM models the latitudinal variations in $\Gamma_{\Omega}$ come from the latitudinal variations in opacity when we discard the small correcting function $\zeta(\theta)$. As a consequence, if the surface opacity were constant (e.g. with Thomson opacity of electrons), $\Gamma_{\Omega}$ would reach unity at all latitudes at the same time when $\omega$ increases! The $\omega$-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 $\omega$-model shows that opacity variations over the stellar surface are unimportant for determining the latitude where $g_{\rm 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_{\rm eff}^{4}/g_{\rm eff}$.

In line with the $\omega$-model and the maximum of $\Gamma_{\Omega}$ at equator, the condition giving the critical angular velocity $\Omega_{c}$ is

[TABLE]

or equivalently,

[TABLE]

These equations show that $\Omega_{c}$ is reduced with increasing Eddington parameter compared to $\Omega_{k}$. Maeder & Meynet (2000) came to the same conclusion, but with a different expression for critical angular velocity, namely $\Omega_{c}\propto\Omega_{k}\sqrt{1-\Gamma_{eq}}$. Because $\Gamma_{eq}\leq 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 $\Gamma_{eq}\ll 1$, we find $\Omega_{c}\simeq\Omega_{k}$, as expected. This is also the solution of Maeder & Meynet (2000) for critical angular velocity in this regime.

2.3.5 Some conclusions from the $\omega$-model

The analysis based on the $\omega$-model underlines three important points:

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

This balance of forces is always first reached at the equator when $\Omega/\Omega_{k}$ increases because the ratio $T_{\rm eff}^{4}/g_{\rm eff}$ at the equator diverges when criticality is reached. 3. 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 $\omega$-model immediately raise the question of the accuracy of this model. This is the next point that we discuss in light of observations and full 2D ESTER models.

2.4 The $\Omega\Gamma$-limit with ESTER 2D models

2.4.1 Interferometric observations (again)

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

2.4.2 Accuracy of the $\omega$-model

The assumptions of the $\omega$-model are that the flux vector is anti-parallel to the effective gravity, the gravitational field is that of a point mass (the Roche model), and the rotation is uniform. The last two of these approximations probably entail the largest errors. They can be appreciated by comparing the flux latitudinal distribution of the $\omega$-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 $\omega$-model, Eq. (11), where $L$, $M$, $\omega$, $r,$ and $g_{\rm eff}$ were taken from the output of ESTER 2D models. For two 2D ESTER ZAMS models of 15 M⊙ and 40 M⊙, we computed the relative difference between the fluxes of the ESTER and $\omega$-model, namely

[TABLE]

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\leavevmode\nobreak\ \%$. For angular velocity ratios lower than $50\leavevmode\nobreak\ \%,$ this difference drops to less than one percent, making the $\omega$-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 $\omega$-model probably never exceeds 10 % for stars with $\omega\leq 0.9$ during the MS. This is illustrated in Fig. 2 (top) with a 15 M⊙ model at mid-MS rotating with $\omega=0.9$. Clearly, the relative difference is reduced compared to the ZAMS model.

2.4.3 The $\Omega\Gamma$-limit with ESTER 2D models

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

In Figure 3 we illustrate the latitudinal variations in $\Gamma_{\Omega}$ for a $15\leavevmode\nobreak\ M_{\odot}$ ESTER model and for a $40\leavevmode\nobreak\ M_{\odot}$ ESTER model both taken at ZAMS, computed for the metallicity $Z=0.02$ and for various values of $\omega$ that we now define as $\omega=\Omega_{eq}/\Omega_{k}$ (see below).

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

To strengthen the effects of radiative acceleration, we considered the case of a $40\leavevmode\nobreak\ M_{\odot}$ ZAMS model. Here we also see (Fig. 3 bottom) that $\omega$ must be as high as $\sim 0.96$ for the $\Omega\Gamma$ limit to be reached222 $\omega$ has to be slightly lower than $0.96$ so that $\Gamma_{\Omega}$ is exactly unity at equator. At $\omega=0.96,$ the star is already supercritical at the equator (grey area in Fig. 3).. 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 $\Omega_{c}$ and $\Omega_{k}$.

From Eq. (28), at ZAMS, we find the equatorial Eddington parameter $\Gamma_{eq}\simeq 0.033$ at criticality for the $15\leavevmode\nobreak\ M_{\odot}$ ESTER model and $\Gamma_{eq}\simeq 0.18$ for the $40\leavevmode\nobreak\ M_{\odot}$ ESTER model. This is surprisingly low for such massive stars. To clarify this result, Fig. 4 shows the latitudinal variations in Eddington parameter $\Gamma$ for both the $15\leavevmode\nobreak\ M_{\odot}$ and $40\leavevmode\nobreak\ M_{\odot}$ ESTER 2D models and for various angular velocity ratios. For the two models, $\Gamma$ decreases with co-latitude when rotation is non-zero. The more rapid the rotation, the lower $\Gamma_{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 $\kappa$ are somewhat unimportant for determining the spatial location of criticality, they are crucial to the value of $\Omega_{c}/\Omega_{k}$.

It might be wondered, however, whether more evolved or more massive stars might have larger $\Gamma_{eq}$ and thus a critical angular velocity that is farther from the Keplerian angular velocity. Fig. 5 shows the evolution of both $\Gamma_{\Omega}(\pi/2)$ and $\Gamma_{eq}$ as a function of the fractional abundance of hydrogen in the convective core $X_{\rm core}/X_{0}$ for a $15\leavevmode\nobreak\ M_{\odot}$ ESTER 2D model initially rotating at $\omega_{i}\equiv\Omega_{eq,i}/\Omega_{k}=0.5$, and without considering any mass loss. For the non-rotating case, evolution tends to increase $\Gamma_{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 $\omega_{i}=0.5$, criticality is reached when $X_{\rm core}/X_{0}\simeq 0.36$. At this time, $\Gamma_{\Omega}$ diverges. While evolution proceeds and $\omega$ grows, the star flattens considerably, causing a significant drop in opacity in the equatorial region. This is clearly shown by the $\Gamma_{eq}$ curve of Fig. 5. After a slight increase at the beginning of evolution, $\Gamma_{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\leavevmode\nobreak\ \mbox{\rm M{}_{\odot}} 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 $\Omega$ as a fraction of the equatorial Keplerian angular velocity $\Omega_{k}$ to appreciate the distance to criticality, bearing in mind that this fraction is slightly smaller than the actual one.

3 Local mass-flux prescription from 1D models

We now address the second question of the paper: the dependence of mass and angular momentum losses on rotation rate. All hot stars have radiation-driven winds that become directly observable in spectral energy distributions and spectral lines as soon as they are above some luminosity threshold. For massive stars of spectral types O, B, and A, this threshold corresponds to $L\sim 10^{4}L_{\odot}$ (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$\alpha$, see Abbott 1979 or Kudritzki & Puls 2000 and references therein). The radiative acceleration has a significant effect on the total gravity, and as shown in the previous section, reduces the critical angular velocity.

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

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

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

[TABLE]

and the radial momentum equation,

[TABLE]

where $\dot{M}$ is the total mass-loss rate of the star, and $\rho$, $v$, $p$ are the density, radial velocity, and gas pressure, respectively. $g_{\rm rad}$ is the radiative acceleration,

[TABLE]

where $g_{\rm rad}^{\rm line}$ is the line-driven acceleration, and $g_{e}=\kappa_{e}F/c$ is the radiative acceleration due to Thomson scattering. $\kappa_{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,

[TABLE]

Using the ideal gas equation of state, we write

[TABLE]

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)

[TABLE]

where $M(t)=kt^{-\alpha}$ is the CAK force multiplier. $\alpha$ and $k$ are the CAK force multiplier parameters (FMPs). $\alpha$ can be interpreted as the ratio of the line force from optically thick lines to the total line-force, which thus decreases with decreasing effective temperature because of the increased iron group lines (e.g. Puls et al., 2000). Moreover, the quantity $k$ is related to the fraction of the total stellar flux, which would be blocked in the photosphere if all lines were optically thick (Puls et al., 2000). $t$ is the electron optical depth parameter and $v_{\rm th}$ is the thermal speed, usually taken as the proton thermal speed $v_{\rm th}\equiv(2k_{\rm B}T_{\rm eff}/m_{\rm 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_{\rm th}\equiv(2k_{\rm B}T_{\rm eff}/m_{\rm Fe})^{1/2}$ for the standard CAK formalism. We then obtain the global mass-loss rate of a non-rotating star, namely

[TABLE]

The global mass loss (without rotation) thus scales as

[TABLE]

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

3.2 Finite disc and ionisation corrections

In the CAK approach, the purely radial streaming radiation leads to an electron optical depth parameter $t$ that only depends on $(dv/dr)^{-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

[TABLE]

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

[TABLE]

where $w=v/v_{\rm th}$ and $w^{\prime}=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}\simeq R$. This assumption may not be verified for rotators close to criticality, for which the fast-wind solution (with $r_{c}\simeq R$) is replaced by the so-called $\Omega$–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 $\Omega$–slow solution was quite modest (factor $\sim 2$ in Curé 2004), we decided to ignore it. We further assumed the velocity to follow a power law like

[TABLE]

where $v_{\infty}$ is the terminal velocity of the wind and 0.7\raisebox{-3.01385pt}{\stackrel{{\scriptstyle<}}{{,\sim,}}}\beta\raisebox{-3.01385pt}{\stackrel{{\scriptstyle<}}{{,\sim,}}}1.3. The corrected force multiplier then simplifies to

[TABLE]

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

[TABLE]

Because $0<\alpha<1$, the finite-disc correction reduces the global mass-loss rate, and for typical values of $\alpha$, namely between 0.4 and 0.7, the CAK mass-loss rate is multiplied by a factor $\sim 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)^{\delta}$ (Abbott, 1982), where $n_{e}$ is the electron density in units of $10^{11}\rm cm^{-3}$ and $W\equiv 0.5(1-\sqrt{1-u^{2}})$ is the radiation dilution factor. $\delta$ is then another FMP. This modification of the line-driven acceleration can be roughly accounted for by replacing $\alpha$ in the power exponents of Eqs. (36) and (42) with $\alpha^{\prime}\equiv\alpha-\delta$ (Puls et al., 1996, 2000). Finally, we obtain the modified local mass-flux in the non-rotating case,

[TABLE]

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 $\kappa_{e}F/c$.

3.3 Parametrisation of the FMPs

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

For $\alpha,$ we took fixed values at $T_{\rm eff}=10\leavevmode\nobreak\ {\rm kK},\ 20\leavevmode\nobreak\ {\rm kK},\ 30\leavevmode\nobreak\ {\rm kK,}\ {\rm and}\ 40\leavevmode\nobreak\ {\rm kK}$ and imposed linear interpolation in between (J. Puls priv. comm.), namely

[TABLE]

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_{\text{eff}}\simeq 25\leavevmode\nobreak\ {\rm kK}$, the mass-loss rate $\dot{M}$ suddenly increases (towards lower $T_{\rm 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_{\text{eff}}=10\leavevmode\nobreak\ {\rm 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 $\dot{M}$ 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_{\rm eff}^{\rm jump}\simeq 20\leavevmode\nobreak\ {\rm kK}$, while Vink et al. (2001) predicted an $\dot{M}$-jump by a factor $\sim 10$ located at ${T_{\rm eff}^{\rm jump}}\simeq 25\leavevmode\nobreak\ {\rm kK}$. Crowther et al. (2006), on the other hand, found a more gradual decrease in terminal velocity $v_{\infty}$ instead (thus a more gradual increase in $\dot{M}$). In addition, a discrepancy of a factor $2$–$3$ also appears when the mass-loss rates of hot OB stars are compared with $T_{\rm eff}>T_{\rm eff}^{\rm jump}$ obtained with the Vink et al. (2001) models and from X-ray, UV, and IR diagnostics (e.g. Najarro et al., 2011; Sundqvist et al., 2011; Bouret et al., 2012; Cohen et al., 2013; Leutenegger et al., 2013; Hervé et al., 2013; Rauw et al., 2015). This could be due to the significant effect of small-scale inhomogeneities in the wind (e.g. Puls et al., 2008; Sundqvist & Owocki, 2013; Puls et al., 2015) and/or to the outdated solar mixture used in the Vink et al. (2001) models, namely $Z_{\odot}\simeq 0.02$ with an Anders & Grevesse (1989) mixture. The more recent solar composition with $Z_{\odot}\simeq 0.014$ of Asplund et al. (2009) could reduce the discrepancy between predicted mass-loss rates and observations (see Section 4.3 for a short discussion of the effects of metallicity on mass loss). Nevertheless, the Vink et al. (2001) models are still widely used in stellar evolution codes, and we also used their recipe to calibrate $k$ to qualitatively predict the impact of radiation-driven winds on the rotational evolution of massive stars.

We thus assumed $\dot{M}=\dot{M}_{\rm 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_{\rm th}$ as the thermal velocity of Fe ions, namely $v_{\rm th}\equiv(2k_{\rm B}T_{\rm eff}/m_{\rm Fe})^{1/2}$. Because we calibrated $k$ using the Vink et al. (2001) mass-loss prescription, each line was considered with its appropriate $v_{\rm th}$.

From now on, the calibrated $k$ is written $k^{\prime}$. We find that $k^{\prime}$ 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^{\prime}$ finally gives us the following semi-empirical function $k^{\prime}(T_{\rm eff})$ at $Z=0.02$, defined on both sides of the effective temperature jump,

[TABLE]

where Vink et al. (2001) defined

[TABLE]

with $\left\langle\rho\right\rangle$ the characteristic wind density at 50 % of the terminal velocity of the wind, given by

[TABLE]

The function $k^{\prime}(T_{\rm eff})$ is shown in Fig. 7. Our values of $k^{\prime}$ assume $v_{\rm th}=(2k_{\rm B}T_{\rm eff}/m_{\rm Fe})^{1/2}$; other assumptions on $v_{\rm th}$ would lead to other values of $k^{\prime}$ to remain compatible with the Vink et al. (2001) mass loss.

4 Effects of rotation on mass and angular momentum loss

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

[TABLE]

This local mass-flux expression is now $\theta$ dependent and thus leads to an anisotropic stellar wind that, at first glance, would favour polar ejection due to the higher polar radiative flux333Our approach implicitly assumes the presence of a weak, polewards directed component of the radiation force. Such a non-radial component is the result of the decreasing radial velocity towards the equator, and is essential for inhibiting a flow that otherwise would be directed towards the equator. Within our approach, however, this component can be neglected when estimating the theta-dependence of $\dot{m}$. For details, see Owocki et al. (1998b), for example. (e.g. Owocki & Gayley, 1997; Owocki et al., 1998b; Petrenz & Puls, 2000; Maeder & Meynet, 2000). We note that Curé (2004) took a different approach and analytically derived an equation for the mass-loss rate that accounts for rotation at the equator. This equation has a $\Omega$–slow solution for rotators close to criticality. We now investigate the surface distribution of $\dot{m}$ from the outputs of ESTER 2D models at ZAMS with $Z=0.02$ and for various $\omega$.

4.1 Latitudinal variations in mass and angular momentum loss

We computed the local mass-flux $\dot{m}(\theta)$ as well as the local angular momentum flux

[TABLE]

with ESTER 2D models and prescription (48).

Fig. 8 shows the local mass-flux $\dot{m}$ and local angular momentum flux $\dot{\ell}$ as a function of colatitude for ESTER 2D models of a $15\leavevmode\nobreak\ M_{\odot}$ star with $Z=0.02$, at ZAMS and for various angular velocity ratios $\omega=\Omega_{eq}/\Omega_{k}$. When the star is far from the $\Omega$-limit, mass loss is favoured in polar regions and thus decreases towards the equator. However, for a sufficiently high angular velocity ratio ($\omega\gtrsim 0.85$ in Fig. 8), there is a colatitude $\theta_{\rm jump}$ where $T_{\text{eff}}(\theta_{\rm jump})=T_{\text{eff}}^{\text{jump}}$ and the local bi-stability limit is crossed. In that case, the mass flux is enhanced between $\theta_{\rm jump}$ and the equator and the star is in a two-wind regime (TWR), otherwise it is in a single-wind regime (SWR). This local bi-stability jump therefore strongly modifies the distribution of the mass flux with colatitude: while mass loss is favoured in polar regions in the SWR, it is far stronger in equatorial regions in the TWR. Similarly, while the angular momentum flux is maximum at some intermediate colatitude in the SWR, it is strongly favoured in the equatorial regions in the TWR. The idea of an enhanced equatorial mass-flux that is due to both gravity darkening and the local bi-stability limit has also been discussed in the past, for instance, by Zickgraf et al. (1986), Zickgraf et al. (1989), Lamers & Pauldrach (1991), Owocki & Gayley (1997), Owocki et al. (1998a), 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 $\omega$ 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 $\omega$.

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

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

[TABLE]

[TABLE]

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

[TABLE]

where $R_{\theta}=\partial R/\partial\theta$ (Rieutord et al., 2016).

The global mass-loss rate $\dot{M}$, the critical angular velocity ratio $\omega_{c}=\Omega_{eq}/\Omega_{c}$ as given by the $\omega$-model (Eq. 28), the ratio of equatorial angular velocity to Keplerian angular velocity $\omega=\Omega_{eq}/\Omega_{k}$ , and the angular momentum loss timescale ${\rm T}_{L}=\mathcal{L}/\dot{\mathcal{L}}$ are reported in Table 1 for a $15\leavevmode\nobreak\ M_{\odot}$ 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 $\omega$, for instance, $\dot{M}(\omega=0.6)/\dot{M}(\omega=0)\simeq 0.87$ (see Fig.9 top). Similar results have been obtained by Müller & Vink (2014).

With increasing $\omega$, the total angular momentum of the star $\mathcal{L}$ increases, and even though the global mass-loss rate $\dot{M}$ decreases in the SWR, the global loss of angular momentum $\dot{\mathcal{L}}$ also increases in this regime. This is simply because $\dot{\mathcal{L}}$ increases for increasing $\omega$. It is even more interesting, however, that the timescale of angular momentum loss ${\rm T}_{\mathcal{L}}=\mathcal{L}/\dot{\mathcal{L}}$ is approximately independent of the degree of criticality $\omega$ in the SWR (see Figure 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, ${\rm T}_{\mathcal{L}}$ rapidly decreases as $\omega$ approaches unity. Both the strong increase in $\dot{M}$ and decrease in ${\rm T}_{\mathcal{L}}$ can be explained by the increasing stellar surface fraction where the effective temperature is lower than $T_{\rm eff}^{\rm jump}$ as $\omega$ increases (see Fig. 8).

That $\dot{M}$ only gradually increases with increasing $\omega$ in the TWR is a result specific to 2D models. In 1D models, the bi-stability jump is accounted for with a stronger global mass-loss rate if the mean effective temperature of the star is lower than $T_{\rm eff}^{\rm jump}\sim 22.5$–$25\leavevmode\nobreak\ {\rm kK,}$ according to Vink et al. (2001). In the present work however, 2D models reach the bi-stability limit if the local effective temperature somewhere on the stellar surface is lower than $T_{\rm eff}^{\rm jump}$. This difference has two consequences. Firstly, 2D models can reach the bi-stability limit even with an average effective temperatures higher than $T_{\rm eff}^{\rm jump}$. In Fig. 10, we illustrate the global mass-loss rate for a variety of angular velocity ratios at ZAMS for a 15 $M_{\odot}$ and a $10\leavevmode\nobreak\ M_{\odot}$ ESTER model at ZAMS and with $Z=0.02$, against the corresponding surface-averaged effective temperature $\overline{T}_{\rm eff}$ of the model. In these models, $\overline{T}_{\rm eff}$ is greater than $T_{\rm eff}^{\rm jump}\simeq 22.8\leavevmode\nobreak\ {\rm kK}$ for all $\omega$. 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_{\rm eff}^{\rm jump}$ monotonically increases with increasing $\omega$. This results in a gradual variation in $\dot{M}$ (and ${\rm T}_{\mathcal{L}}$) with $\overline{T}_{\rm eff}$ in the TWR (see Fig. 10). Hence, in rotating stars the bi-stability jump does not induce a discontinuity of the global mass-loss rate (but it induces a discontinuity of the local mass-flux, see Fig. 8) as the $\omega$ parameter increases (and therefore as the mean effective temperature decreases). The discontinuity occurs only on the derivative of the function $\dot{M}(T_{\rm 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 line-broadening effects such as the so-called macro-turbulence (e.g. Simón-Díaz & Herrero, 2007).

4.3 Metallicity effect

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

5 Discussion and conclusions

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

To the first question, we answer that the critical angular velocity is very close to the classical Keplerian angular velocity at the equator, at least for stars with masses lower than 40\leavevmode\nobreak\ \mbox{\rm M{}_{\odot}} (and for $Z=0.02$). The role of radiative acceleration turns out to be rather limited because of the combination of a reduced opacity and reduced flux at the equator. The reduction of the flux, the so-called gravity darkening, is less than was predicted by the von Zeipel model. This latter point is the main difference between this study and the pioneering investigations of Maeder (1999) and Maeder & Meynet (2000). ESTER 2D models indeed show that the flux is almost anti-parallel to gravity in the stellar radiative envelope (Espinosa Lara & Rieutord, 2011). To a very good approximation, we can therefore write $\@vec{F}=f(r,\theta)\@vec{g}_{\rm eff}$, which is the base of the $\omega$-model (Espinosa Lara & Rieutord, 2011; Rieutord et al., 2016). We showed that the $\omega$-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,\theta)$ 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,\theta)$ 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 $\@vec{g}_{\rm tot}=\@vec{g}_{\rm rad}+\@vec{g}_{\rm 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 $\@vec{g}_{\rm tot}$ to vanish somewhere at the surface, is therefore always reached before the equatorial rotation reaches Keplerian rotation. This point has been made by Maeder & Meynet (2000). However, 2D models hold that this difference is tiny. Considering a massive star of 40 $M_{\odot}$, we therefore find that criticality, $\@vec{g}_{\rm tot}=\@vec{0}$ at the equator, is reached at $\omega\sim 0.96$ and even at 0.997 for a 15 $M_{\odot}$ star. This tiny difference can be understood because gravity darkening in the $\omega$-model is weaker than in the von Zeipel model and because the effect of rotation on opacity leads to a strong decrease in standard Eddington parameter towards the equator. To return to the debate between Langer (1997, 1998) and Glatzel (1998), our results support the latter concerning the influence of the Eddington limit on the value of critical rotation: this influence is quite small and never exceeds 4% as far as we could test (i.e. M$\leq$ 40 $M_{\odot}$, $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 $\dot{M}$ in the range $T_{\rm eff}\in[10,50]{\rm kK}$. This prescription led to a discontinuity in the mass flux when $T_{\rm 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\leavevmode\nobreak\ \mbox{\rm M{}{\odot}}), but their realism and their reliability decrease with increasing mass. Beyond 40\leavevmode\nobreak\ \mbox{\rm M{}{\odot}}, 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 $\dot{M}(T_{\rm 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_{\rm eff}$, implying that any discontinuity in the $\dot{M}(T_{\rm eff})$ relation is smoothed out, unless we can reproduce the observed star with a complete 2D wind+star model. In a follow-up paper (Gagnier et al., 2019), we apply these results to study the evolution of rotation in early-type fast-rotating stars and address the question, among others, how a wind can prevent a massive star from reaching the critical rotation.

Acknowledgements.

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

Appendix A Opacity dependence on effective temperature

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

[TABLE]

where $\rho_{s}(\theta)$ is the local surface density. This approximation seems to be rather well verified at the surface of rapidly rotating ESTER 2D models. We recall that the pressure at the surface is

[TABLE]

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

[TABLE]

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

Appendix B Short presentation of ESTER models

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

[TABLE]

where $\phi$ is is the gravitational potential; the continuity equation,

[TABLE]

the momentum equation,

[TABLE]

where $\@vec{F}_{\rm visc}$ is the viscous force; and the heat balance equation,

[TABLE]

and

[TABLE]

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

These equations are completed by boundary conditions that require that (i) the gravitational potential $\phi$ vanishes at infinity, (ii) velocity fields meet stress-free conditions at the stellar surface, (iii) that the surface radiates like a local black body, and (iv) the surface is defined by the place where the pressure $P$ equals $g_{\rm eff}/\kappa$. Usual notations have been used: $G$ is the gravitational constant, $\@vec{v}$ the velocity field, $s$ the entropy, and $\varepsilon_{*}$ the energy produced by nuclear reactions per unit mass.

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