Issue 
A&A
Volume 646, February 2021



Article Number  A180  
Number of page(s)  12  
Section  Stellar atmospheres  
DOI  https://doi.org/10.1051/00046361/202039224  
Published online  26 February 2021 
Analytic, dustindependent massloss rates for red supergiant winds initiated by turbulent pressure
^{1}
Institute of Astronomy,
KU Leuven, Celestijnenlaan 200D,
3001
Leuven,
Belgium
email: dylan.kee@kuleuven.be
^{2}
Anton Pannekoek Institute, University of Amsterdam,
Science Park 904,
1098XH
Amsterdam, The Netherlands
Received:
20
August
2020
Accepted:
19
December
2020
Context. Red supergiants are observed to undergo vigorous mass loss. However, to date, no theoretical model has succeeded in explaining the origins of these objects’ winds. This strongly limits our understanding of red supergiant evolution and Type IIP and IIL supernova progenitor properties.
Aims. We examine the role that vigorous atmospheric turbulence may play in initiating and determining the massloss rates of red supergiant stars.
Methods. We analytically and numerically solve the equations of conservation of mass and momentum, which we later couple to an atmospheric temperature structure, to obtain theoretically motivated massloss rates. We then compare these to stateoftheart empirical massloss rate scaling formulae as well as observationally inferred massloss rates of red supergiants.
Results. We find that the pressure due to the characteristic turbulent velocities inferred for red supergiants is sufficient to explain the massloss rates of these objects in the absence of the normally employed opacity from circumstellar dust. Motivated by this initial success, we provide a first theoretical and fully analytic massloss rate prescription for red supergiants. We conclude by highlighting some intriguing possible implications of these rates for future studies of stellar evolution, especially in light of the lack of a direct dependence on metallicity.
Key words: stars: massloss / stars: winds, outflows / stars: massive / supergiants / turbulence
© ESO 2021
1 Introduction
The lack of a satisfactory theory explaining the strong, >10^{−7} M_{⊙} yr^{−1}, mass loss for evolved massive stars on the red supergiant branch has been a longstanding problem in our understanding of these objects (see Levesque 2017, for a recent discussion of this and other current puzzles in studying red supergiants). While this considerable gap in our knowledge has been patched over somewhat by empirical rates and scaling formulae (see, e.g., Mauron & Josselin 2011, for a review), the overall disagreement between any given two of these leaves those attempting to model red supergiants directly, or attempting to model stellar evolution including or approaching this phase, with a rather untenable problem. Namely, one must somehow distinguish and select between rates that may differ from each other by as much as four orders of magnitude in some parts of the parameter space without a general understanding of the underlying process that determines this mass loss.
In order to understand why theories generally struggle so much in this region, it is important to highlight some general features of previous modeling attempts. In analogy to the lowermass asymptotic giant branch (AGB) stars, it has generally been assumed that the primary driving mechanism of red supergiant winds is radiation pressure on dust grains (see, e.g., Höfner & Olofsson 2018, for a review). For AGB stars, pulsations provide an atmospheric piston levitating material off the stellar surface and up to a region where the temperature has dropped far enough to allow dust to condense (e.g., Bladh & Höfner 2012; Bladh et al. 2013). However, as dust nucleation is also a densitydependent process, it is essential not only that material reaches this region but that enough material is levitated to make the dust formation efficient. Due to both the lower pulsational amplitudes inferred for red supergiants (e.g., Wood et al. 1983) and their higher effective temperatures relative to AGB stars, applying similar models to red supergiants thus requires gas to reach a greater height in comparison to the stellar radius, while less material is expected to be levitated in the first place. This implies a strongly decreased efficiency of dust condensation. Modeling attempts have been generally unsuccessful in generating the atmospheric extensions of red supergiants necessary to put enough material at the dust sublimation front to recover the observed massloss rates of red supergiants (e.g., ArroyoTorres et al. 2015).
An alternative suggestion has been that pulsational motions might be accompanied or replaced by significant atmospheric turbulence (e.g., Gustafsson & Plez 1992; Josselin & Plez 2007), and that this turbulence might be seeded by the vigorous convection expected in the atmospheres of red supergiants (see Freytag et al. 2012, for a theoretical treatment of convection in red supergiant envelopes). At the moment, these convection simulations produce only modestly extended atmospheres compared to what is inferred from observations (e.g., ArroyoTorres et al. 2015). Nevertheless, observations of red supergiants do suggest that the outer layers of these stars are indeed very turbulent (e.g., Josselin & Plez 2007; Ohnaka et al. 2017).
Therefore, inspired by the work of Gustafsson & Plez (1992) and Josselin & Plez (2007), we here undertake the derivation of an analytic, theoretical massloss rate that focuses on these observed turbulent velocities present in the atmospheres of red supergiants. Doing so, we find that the inferred turbulent motions are, alone, sufficient to explain the massloss rates of red supergiants even in absence of any dust opacity. This model then can be seen as a turbulent pressure driven extension of the classical thermal pressure driven Parker wind (Parker 1958). In contrast to the solar wind and other applications of this theory (for instance the “warm wind” model, Hearn 1975; Rogerson & Lamers 1975), however, the total sound and turbulent speed can be much lower due to the lower surface gravity of red supergiants. In Sect. 2 we lay out the basic formalism for this model in an isothermal setup, under which the massloss rate is fully analytic. Byimposing a temperature structure, we then extend the model in Sect. 3 by iteratively solving the relevent equations to converge to nonisothermal massloss rates. We further fit over the difference between theisothermal and nonisothermal massloss rates to provide a theoretically motivated mass loss rate as an analytic function of stellar parameters and turbulent velocity. In Sect. 4 we review some of the currently used empirical massloss rate prescriptions for red supergiants and compare our theoretical rate to these, as well as to the sample of observed red supergiants from Josselin & Plez (2007) and the observations of Antares by Ohnaka et al. (2017). Finally, we discuss some future directions for this model in Sect. 5.
2 Analytic massloss rate from levitation of an isothermal atmosphere
2.1 Derivation of the massloss rate
To begin, we first consider the relevant equations that must be satisfied at all points in an isothermal flow, namely conservation of mass and momentum. For a purely radial, spherically symmetric outflow in a steady state these are $$\begin{array}{l}v\frac{\partial \rho}{\partial r}+\rho \frac{\partial v}{\partial r}+\frac{2\rho v}{r}=0\\ v\frac{\partial v}{\partial r}=\frac{1}{\rho}\frac{\partial P}{\partial r}\frac{G{M}_{\ast}}{{r}^{2}},\end{array}$$
expressed in terms of radial velocity v, density ρ, and total pressure P as functions of spherical radius r, as well as stellar mass M_{*} and gravitational constant G.
As is done in standard static model atmosphere computations (see, e.g., Gustafsson et al. 2008; their Sect. 3) pressure is expressed as the sum of thermal pressure P_{therm} = $\rho {c}_{\text{s}}^{2}$, turbulent pressure ${P}_{\text{turb}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\rho {v}_{\text{turb}}^{2}$, and radiation pressure P_{rad} where the characteristic velocities of the first two are respectively the sound speed c_{s} and the turbulent speed v_{turb}. As concerns the radial component of radiation pressure gradient, we associate this with the acceleration g_{rad} such that $$\frac{1}{\rho}\frac{\partial {P}_{\text{rad}}}{\partial r}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{g}_{\text{rad}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{\kappa {F}_{r}}{c}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{\kappa {L}_{\ast}}{4\pi {r}^{2}c}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\Gamma \frac{G{M}_{\ast}}{{r}^{2}},$$(3)
with c the speed of light, κ the flux weighted mean opacity [cm^{2} g^{−1}], and F_{r} the radial component of the flux, related to the stellar luminosity L_{*} as F_{r} = L_{*}∕(4πr^{2}). The final equality further introduces the Eddington factor Γ ≡ κL_{*}∕(4πGM_{*}c). Explicitly replacing pressure with this combination and substituting Eq. (1) into Eq. (2) to eliminate density gradient terms, we find $$v\left(1\frac{{c}_{\text{s}}^{2}+{v}_{\text{turb}}^{2}}{{v}^{2}}\right)\frac{\partial v}{\partial r}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{2\left({c}_{\text{s}}^{2}+{v}_{\text{turb}}^{2}\right)}{r}\frac{G{M}_{\ast}\left(1\Gamma \right)}{{r}^{2}},$$(4)
where we have assumed that v_{turb} is constant. We discuss this approximation further below.
We note that Eq. (4) is the same equation as is solved for an isothermal, pressure driven Parker wind (Parker 1958), only now with a modified “effective sound speed”, ${c}_{\text{s},\text{eff}}\equiv \sqrt{{c}_{\text{s}}^{2}+{v}_{\text{turb}}^{2}}$. At the location where the flow velocity reaches this effective sound speed the left side of Eq. (4) goes to zero, and therefore the right side must also go to zero. Solving for this criteria yields that the location of this critical point in the wind outflow is given by a modified Parker radius, $${R}_{\text{p},\text{mod}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{G\text{\hspace{0.17em}}{M}_{\ast}\text{\hspace{0.17em}}(1\Gamma )}{2\text{\hspace{0.17em}}({c}_{\text{s}}^{2}+{v}_{\text{turb}}^{2})}.$$(5)
This then suggests that the problem of wind acceleration here can be envisioned to consist of two parts. First, a lowspeed near “levitation” of material to this modified Parker radius at which v = c_{s,eff}, and, second, a subsequent acceleration of the material to infinity.
Below the modified Parker radius the contribution of the advection term (v ∂v∕∂r) in Eq. (2) is almost negligible so that the equation reduces to the standard equation for hydrostatic equilibrium, $$\frac{1}{\rho}\frac{\partial \rho}{\partial r}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{G\text{\hspace{0.17em}}{M}_{\ast}\text{\hspace{0.17em}}(1\Gamma )}{({c}_{\text{s}}^{2}+{v}_{\text{turb}}^{2})\text{\hspace{0.17em}}{r}^{2}}.$$(6)
Typically, the Eddington factor and, as alluded to above, the turbulent speed may be expected to vary in radius. However, under the assumption that this variation is only mild, we may replace the exact expression Γ(r) with Γ and, as already done in Eq. (4), v_{turb}(r) with v_{turb}. These representative constant values of Γ and v_{turb} then in general would designate some appropriate spatial means (see also the discussion in Gustafsson et al. 2008). Then Eq. (6) can be integrated analytically from the stellar radius R_{*} to an arbitrary radius r ≤ R_{p,mod} using the boundary value definition ρ_{*}≡ ρ(R_{*}) to give $$\rho (r)\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{\rho}_{\ast}\mathrm{exp}\left[\frac{{R}_{\ast}}{H}\left(1\frac{{R}_{\ast}}{r}\right)\right],$$(7)
where we have made the substitution $$\frac{H}{{R}_{\ast}}\equiv \frac{{R}_{\ast}\text{\hspace{0.17em}}({c}_{\text{s}}^{2}+{v}_{\text{turb}}^{2})}{G\text{\hspace{0.17em}}{M}_{\ast}\text{\hspace{0.17em}}(1\Gamma )}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}2\text{\hspace{0.17em}}\frac{{c}_{\text{s}}^{2}+{v}_{\text{turb}}^{2}}{{v}_{\text{esc},\text{eff}}^{2}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{1}{2}\text{\hspace{0.17em}}\frac{{R}_{\ast}}{{R}_{\text{p},\text{mod}}}.$$(8)
Here ${v}_{\text{esc},\text{eff}}\equiv \sqrt{G{M}_{\ast}(1\Gamma )/{R}_{\ast}}$ denotes the escape speed from R_{*} with $\sqrt{1\Gamma}$ accounting for the reduction in effective gravity due to radiative acceleration.
In order to determine the value of ρ_{*}, we make use of the definition of the stellar radius in optical depth space, τ(R_{*}) ≡ 2∕3, from a spherical extension of the classical planar gray model atmosphere (Lucy 1971; and see further Sect. 3). Within this model, one defines a “spherically modified” optical depth scale as $$\tau (r)\equiv {\displaystyle {\int}_{r}^{\infty}\kappa}\text{\hspace{0.17em}}\rho \text{\hspace{0.17em}}{\left(\frac{{R}_{\ast}}{r}\right)}^{2}\text{d}r,$$(9)
and, as such, $$\kappa {\displaystyle {\int}_{{R}_{\ast}}^{\infty}\rho}(r){\left(\frac{{R}_{\ast}}{r}\right)}^{2}\text{d}r\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{2}{3},$$(10)
where we have again used the above approximation that Γ, and therefore κ, is independent of radius. Inserting the hydrostatic density profile from Eq. (7) into Eq. (9) yields an analytic value of ρ_{*}, $${\rho}_{\ast}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{2}{3\text{\hspace{0.17em}}\kappa \text{\hspace{0.17em}}H}{\left[1\mathrm{exp}\left(\frac{{R}_{\ast}}{H}\right)\right]}^{1}.$$(11)
Before combining the above terms to derive a turbulent pressure driven massloss rate, it is important to examine a few of the assumptions that went into the derivation. First, we consider the approximation that the advection term is negligible for r < R_{p,mod}. If we relax this approximation, we find that integrating Eq. (2) from R_{*} to an arbitrary radius r gives $$\mathrm{log}\left(\frac{\rho (r)}{{\rho}_{\ast}}\right)\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{{R}_{\ast}}{H}\left(1\frac{{R}_{\ast}}{r}\right)\frac{1}{2\text{\hspace{0.17em}}({c}_{\text{s}}^{2}+{v}_{\text{turb}}^{2})}{\displaystyle {\int}_{{R}_{\ast}}^{r}\frac{\partial {v}^{2}}{\partial r}}\text{d}r.$$(12)
While in general solving this equation at an arbitrary radius requires numerical integration, the integration from R_{*} to R_{p,mod} remains analytic as $v{({R}_{\ast})}^{2}/({c}_{\text{s}}^{2}+{v}_{\text{turb}}^{2})\approx 0$ and $v{({R}_{\text{p},\text{mod}})}^{2}/({c}_{\text{s}}^{2}+{v}_{\text{turb}}^{2})\equiv 1$. Comparing ρ(R_{p,mod}) computed with and without the advection term shows that including this term reduces the density at the modified Parker radius by a constant factor exp(−1∕2). As this reduction is analytic and constant, we therefore include it in all calculations for the remainder of this section.
Related to this change in density, it is also important to examine whether the optical depth, and as such the base density ρ_{*}, is significantly impacted by the inclusion of the advection term or by a radial profile in the opacity. To test this we plot τ(r∕R_{*}) computed from the hydrostatic density structure with a constant opacity in Fig. 1, assuming that H∕R_{*} = 0.07, a characteristic value for a red supergiant atmosphere. We note here, from examining Eqs. (7), (9), and (11), that τ(r∕R_{*}) only depends on this chosen H∕R_{*} as the definition of τ(1) = 2∕3 results in κ in Eqs. (9) and (11) cancelling. As expected, almost all of the optical depth is accumulated within only a few scale heights of the stellar surface, here within about half a stellar radius of R_{*}. This is well away from the modified Parker radius, here ~7R_{*}, and therefore the optical depth used to define the stellar radius and ρ_{*} is all accumulated in a region where the advection term is indeed negligible. Similarly, this means that κ and $({v}_{\text{turb}}^{2}+{c}_{\text{s}}^{2})$ only need to be constant over this same small spatial extent in order for ρ_{*} to be unaffected by any variations they may have. Section 3 will revisit the implications of allowing c_{s} to take on a realistic radial profile.
Finally, to turn the terms derived thus far into a massloss rate, we note that Eq. (1) is equivalent to the constraint ρvr^{2} = const.; the total massloss rate is this conserved quantity integrated over all solid angle, Ṁ ≡4πρvr^{2}. Further, as this quantity is independent of radius, we can simply take its value at the modified Parker radius $$\dot{M}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}4\text{\hspace{0.17em}}\pi \text{\hspace{0.17em}}\rho ({R}_{\text{p},\text{mod}})\text{\hspace{0.17em}}\sqrt{{c}_{\text{s}}^{2}+{v}_{\text{turb}}^{2}}\text{\hspace{0.17em}}{R}_{\text{p},\text{mod}}^{2},$$(13)
as the total massloss rate of the star. Combining Eqs. (7), (8), and (11) and accounting for the additional factor exp(−1∕2) from advection gives ρ(R_{p,mod}) to be $$\rho ({R}_{\text{p},\text{mod}})\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{4}{3}\frac{{R}_{\text{p},\text{mod}}}{\kappa {R}_{\ast}^{2}}\frac{\mathrm{exp}\left[\frac{2{R}_{\text{p},\text{mod}}}{{R}_{\ast}}+\frac{3}{2}\right]}{1\mathrm{exp}\left[\frac{2{R}_{\text{p},\text{mod}}}{{R}_{\ast}}\right]}.$$(14)
All quantities in Eqs. (13) and (14) are analytically known within this formalism, such that these expressions together offer a fully analytic massloss rate.
Fig. 1 Optical depth calculated inward to r. Note that most of the optical depth is accumulated within a stellar radius of the star. 
2.2 Scaling of massloss rate with key parameters of the model
In understanding the regimes of stellar parameters where such a massloss rate model is and is not viable, it is important to note that H∕R_{*} = R_{*}∕(2R_{p,mod}) as shown in Eq. (8). This means that the same scale which controls the exponential stratification of density also controls where the critical point lies. Therefore, H∕R_{*} must be relativelylarge compared to, for instance, the Sun in order to prevent the exponential density stratification of the atmosphere from driving Ṁ to zero. At the same time, R_{p,mod} and by extension R_{*} must be large (in absolute units) simply due to the ${R}_{\text{p},\text{mod}}^{2}$ dependence of Ṁ.
The massloss rate derived in Sect. 2.1 is a function of only a few key parameters of the physical setup, namely stellar mass M_{*}, stellar radius R_{*}, sound speed c_{s}, turbulent velocity v_{turb}, and opacity κ. Specifically, $$\dot{M}\propto \frac{1}{\kappa}{\left(\frac{{R}_{\ast}}{H}\right)}^{3}\mathrm{exp}\left(\frac{{R}_{\ast}}{H}\right)\text{\hspace{0.17em}}{R}_{\ast},$$(15)
where H∕R_{*} itself scales as $$\frac{H}{{R}_{\ast}}\propto \frac{{R}_{\ast}\text{\hspace{0.17em}}\left({c}_{\text{s}}^{2}+{v}_{\text{turb}}^{2}\right)}{{M}_{\ast}\text{\hspace{0.17em}}\left(1\Gamma \right)}.$$(16)
We can trade the sound speed dependence for a stellar parameter by assuming for the following discussion that the wind has ${c}_{\text{s}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\sqrt{{k}_{\text{B}}\text{\hspace{0.17em}}{T}_{\text{eff}}/{m}_{\text{H}}}$ with Boltzmann constant k_{B}, mean molecular mass equal to the mass of a hydrogen atom m_{H}, and temperature equal to the stellar effective temperature ${T}_{\text{eff}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{(L/(4\text{\hspace{0.17em}}\pi \text{\hspace{0.17em}}{R}_{\ast}^{2}\text{\hspace{0.17em}}\sigma ))}^{1/4}$, where σ is the StefanBoltzmann constant. To investigate the behavior of this analytic expression, we begin by calculating the massloss from a star with parameters selected to be consistent with a “typical” RSG, namely M_{*} = 10 M_{⊙}, R_{*} = 500 R_{⊙}, T_{eff} = 3500 K (c_{s} = 5.4 km s^{−1}), v_{turb} = 15 km s^{−1}, and κ = 0.01 cm^{2} g^{−1}. Note that here opacity is estimated from the OPAL opacity tables (Iglesias & Rogers 1996) as a lower limit to the total opacity as will be discussed further below. For these parameters, Ṁ = 8.5 × 10^{−7} M_{⊙} yr^{−1}, in the rangeof plausible massloss rates for a red supergiant (see, e.g., Josselin & Plez 2007).
For the remaining discussion in this section, we vary individual model parameters while holding all others to the values provided immediately above to show how Ṁ reacts to such variations. The scaling of this massloss rate with stellar parameters is relatively straightforward as shown by Fig. 2. Specifically, increasing the stellar radius or effective temperature over a reasonable range for red supergiants increases the massloss rate because increasing either of these increases H∕R_{*}. This goes into the exponential term in Eq. (15) which increases faster than the powerlaw terms fall off. Conversely, increasing M_{*} has the opposite effect on H∕R_{*} and therefore increasing stellar mass decreases massloss rate. Note that T_{eff} only appears through the sound speed, and thus only in the combination ${c}_{\text{s}}^{2}+{v}_{\text{turb}}^{2}$. Therefore, it is unsurprising that varying this parameter has a much weaker effect than either stellar mass or radius, a fact we return to later in the discussion.
Next, weexamine the scaling of massloss rate with v_{turb}, shown in Fig. 3. In order to emphasize the importance of v_{turb} as the driving mechanism for this model, we vary v_{turb} from a minimum microturbulence for red supergiants (~3 km s^{−1}) up through the characteristic turbulent velocities inferred by Josselin & Plez (2007). As was the case with the variation of stellar parameters, for much of the range of v_{turb} the dominant effect is on the scale height. As $H/R\propto {v}_{\text{turb}}^{2}$ for v_{turb} ≫ c_{s}, this dependence is quite steep. However, note here that we allow a much larger range of variation in v_{turb} than in the other parameters, which means that we see an additional regime that did not appear in the variation of stellar parameters, namely the asymptotic massloss rate reached at high v_{turb}. Here exp (−R_{*}∕H) begins to vary less rapidly with the increase in H∕R_{*} and the decrease of R_{p,mod} and ρ_{*} with increasing H∕R_{*} takes over. The scaling is truncated at the point where increasing v_{turb} means R_{p,mod} < R_{*}, as the model is no longer meaningful in this regime. If v_{turb} were allowed to vary to arbitrarily low values, a lower asymptote would also appear corresponding to the negligibly low mass loss that could be driven byc_{s} alone.
Finally, we examine the scaling of Ṁ with opacity, shown in Fig. 4. While the first calculation of Ṁ for characteristic red supergiant parameters above assumed the Rosseland mean opacity from OPAL, the acceleration of the wind means that a purely static opacity model may not be the most accurate one (see also Gustafsson & Plez 1992). In fact, the Doppler shifting of spectral lines from atomic and molecular species will increase the characteristic opacity scale of the wind, as discussed in Appendix A. In practice, the computation presented in Appendix A for κ and Γ from these Doppler shifted spectral lines becomes computationally expensive to include in a model such as the one we discuss in the remainder of the paper. Specifically, millions of spectral lines are involved (see especially the tabulations of infrared spectral lines of water, e.g., Barber et al. 2006), with the opacity of each line becoming positiondependent due both to the velocity gradient and the radial temperature variation. In studies of hot, OB star winds, this is circumvented by introducing a distribution function that allows one to analytically approximate the cumulative effects of all these spectral lines. As a first step, we have taken the brute force approach of computing the full sum but only for fixed hydrodynamic structures, thus allowing us to approximate the net effect of this additional radiative driving. By computing the force from infrared spectral line lists of CO (Goorvitch 1994), TiO (McKemmish et al. 2019), and H_{2} O (Barber et al. 2006) as a sample, we find that the net increase in fluxweighted mean opacity may be an order of magnitude or more over the basal OPAL opacity. However, this preliminary study deserves further attention as it is based on analyses of fixed hydrodynamic structures, while line acceleration is a sensitive function of the velocity field. Further, we have only included three molecules here, and also used a simplified radiative transfer treatment based on the socalled Sobolev approximation (see Appendix A), implying the real effect could be a bit different than the 1 dex mentioned above.
Here we take a more general view of the effect of opacity variation, simply allowing Γ to increase from its basal value as implied by the Rosseland mean opacity, Γ_{Ros}. As was the case for v_{turb}, we truncate the plot when R_{p,mod} = R_{*} as this is again the limit at which the model breaks down. To emphasize the complex behavior induced by varying Γ and thus κ, we show the scaling of Ṁ for three models: v_{turb} = 10 km s^{−1} and R_{*} = 500 R_{⊙} (left panel of Fig. 4), v_{turb} = 15 km s^{−1} and R_{*} = 500 R_{⊙} (center panel of Fig. 4), v_{turb} = 15 km s^{−1} and R_{*} = 1000 R_{⊙} (right panel of Fig. 4). Note here that ${\Gamma}_{\text{Ros}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{\kappa}_{\text{R}}\text{\hspace{0.17em}}{L}_{\ast}/(4\text{\hspace{0.17em}}\pi \text{\hspace{0.17em}}G\text{\hspace{0.17em}}{M}_{\ast}\text{\hspace{0.17em}}c)\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{\kappa}_{\text{R}}\text{\hspace{0.17em}}{R}_{\ast}^{2}\text{\hspace{0.17em}}\sigma \text{\hspace{0.17em}}{T}_{\text{eff}}{}^{4}/(G\text{\hspace{0.17em}}{M}_{\ast}\text{\hspace{0.17em}}c)$. All parameters are the same for the three cases except R_{*}, such that this minimal point for Γ is lower for the two models with R_{*} = 500 R_{⊙} (Γ_{Ros} = 2.6 × 10^{−3}) than the one with R_{*} = 1000 R_{⊙} (Γ_{Ros} = 1.0 × 10^{−2}).
In general, as increasing Γ means increasing the radiation force available to levitate the atmosphere, one might expect that Ṁ would monotonically increase with Γ. However, it is important to recall that we choose our unique solution for the massloss rate by requiring τ(R_{*}) = 2∕3 through Eq. (10), which results in an inverse dependence of the stellar photosphere density on opacity in Eq. (11). When Γ is small and the wind is predominantly driven by turbulent pressure, this inverse dependence wins out and increasing Γ and therefore κ implies a reduction of the density of the wind and therefore also of the massloss rate to keep R_{*} fixed.
As Γ approaches unity, radiation then contributes more meaningfully to the total force budget of the wind and the expected increase manifests itself. In the case that the turbulent pressure gradient is not already driving a strong wind mass loss (e.g. v_{turb} and R_{*} small as in the left panel of Fig. 4) then radiative acceleration is able to drive a much stronger mass loss near Γ unity than pressure could at low Γ. In the case where turbulent pressure already drives a strong mass loss (i.e., v_{turb} and R_{*} large as in the right panel of Fig. 4) the net effect is only to cancel out the reduction in Ṁ from the increased opacity. However, even this increase (or flattening) is not able to hold out all the way to Γ unity. Eventually, as was the case with increasing v_{turb}, R_{p,mod} is driven to R_{*} and Ṁ drops off. Again, this point merely denotes where such a turbulent pressure driven wind model breaks down.
As mentioned in the introduction, it is intriguing that this model is able to recover reasonable massloss rates for red supergiants without appealing to dust opacity. In fact, as Γ_{Ros} is small in the standard case we consider, this model does not appeal to significant radiative driving from any source. One of the theories present in the literature is that red supergiant winds behave analagously to the winds of AGB stars with some atmosphericextension allowing for the condensation of dust, the opacity of which is the main massloss rate driver (see e.g. Höfner & Olofsson 2018, for a review). Instead, we here find that it is plausible that the extreme extension of their atmospheres by turbulent pressure alone could be the dominant part determining the red supergiant massloss process.
However, it is important to note here that opacity may still play an important role in altering the structure of the winds of red supergiants. While we have here assumed a constant opacity as a function of distance from the star in order to examine the scaling of Ṁ with Γ, this is not expected to be the case. As the wind cools away from the star, additional molecules and dust will form, thereby enhancing the wind opacity. Even the continuum opacity is likely to shift as the Hydrogen anion is no longer the dominant continuum opacity source. As discussed in the prior subsection, the effect of this enhanced opacity on at least the optical depth scale, and as such the density scale used in this model, may not be substantial as long as the opacity enhancement occurs beyond the first stellar radius or so. This is because, as shown by Fig. 1, the majority of optical depth is accumulated in a region quite near to the star making κ in Eq. (11) effectively an average over this near star region. Enhancing the opacity further out would, however, contribute to extending the scale height of the atmosphere, and even potentially allow the wind to switch from being turbulence driven to radiation driven at some location. Both these effects could increase Ṁ somewhat over what is presented here. In light of these open items related to wind opacity, as we discuss further in Sect. 5, future work should reexamine what role depth dependent dust, molecular, and continuum opacity may play, especially in setting the terminal velocity of the stellar wind as this is illdefined under a Parkerlike pressure driven wind (see also the discussion in Sect. 3). For the remainder of this paper, however, given the highly promising results throughout, we proceed with turbulent pressure supplemented only by the marginal additional contribution of a constant Γ_{Ros}.
Fig. 2 Dependence of the analytic massloss rate on stellar mass (left), stellar radius (center), and temperature (right). 
Fig. 3 Dependence of the analytic massloss rate on turbulent velocity. 
Fig. 4 Dependence of the analytic massloss rate on Γ for M_{*} = 10 M_{⊙} and T_{eff} = 3500 K. Left panel: uses R_{*} = 500 R_{⊙} and v_{turb} = 10 km s^{−1}, middle panel: uses R_{*} = 500 R_{⊙} and v_{turb} = 15 km s^{−1}, and right panel: uses R_{*} = 1000 R_{⊙} and v_{turb} = 15 km s^{−1}. 
3 Numerically determined massloss rate from steadystate, nonisothermal, dynamical winds
3.1 Iteration procedure to find a unique massloss rate
While the simplified, isothermal treatment has the benefit of being purely analytic, we can also treat the nonisothermal wind structure in a numerical approach. To do so, we use the spherically symmetric generalization to the plane parallel gray atmosphere introduced by Lucy (1971) $$\begin{array}{ll}{T}^{4}\hfill & \text{\hspace{0.17em}}=\text{\hspace{0.17em}}{T}_{\text{eff}}{}^{4}\left(W+\frac{3}{4}{\displaystyle {\int}_{{R}_{\ast}}^{r}\kappa}\text{\hspace{0.17em}}\rho \text{\hspace{0.17em}}{\left(\frac{{R}_{\ast}}{r}\right)}^{2}\text{d}r\right)\hfill \\ \hfill & \text{\hspace{0.17em}}=\text{\hspace{0.17em}}{T}_{\text{eff}}{}^{4}\left(W+\frac{3}{4}\tau \right),\hfill \end{array}$$
where $W\equiv (1\sqrt{1{({R}_{\ast}/r)}^{2}}\text{\hspace{0.17em}})/2$ is the geometric dilution factor. The second equality again uses the spherically modified optical depth as given by Eq. (9). For low optical depth material, and for distances away from the stellar surface this temperature structure behaves like $T\propto 1/\sqrt{r}$ as is generally assumed for the dust free regions near red supergiants (e.g. Decin et al. 2006). Accounting for this temperaturestructure in Eqs. (1) and (2) while still maintaining the assumption that v_{turb} is constant leads to the system of nonisothermal equations $$\begin{array}{l}v\text{\hspace{0.17em}}\frac{\partial \rho}{\partial r}+\rho \text{\hspace{0.17em}}\frac{\partial v}{\partial r}+\frac{2\text{\hspace{0.17em}}\rho \text{\hspace{0.17em}}v}{r}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}0\\ v\text{\hspace{0.17em}}\frac{\partial v}{\partial r}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{{c}_{\text{s}}^{2}+{v}_{\text{turb}}^{2}}{\rho}\text{\hspace{0.17em}}\frac{\partial \rho}{\partial r}\frac{G\text{\hspace{0.17em}}{M}_{\ast}\text{\hspace{0.17em}}(1\Gamma )}{{r}^{2}}\frac{{k}_{\text{B}}}{{m}_{\text{H}}}\text{\hspace{0.17em}}\frac{\partial T}{\partial r}\\ \frac{\partial {T}^{4}}{\partial r}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{T}_{\text{eff}}{}^{4}\left(\frac{\partial W}{\partial r}+\frac{3}{4}\kappa \rho {\left(\frac{{R}_{\ast}}{r}\right)}^{2}\right).\end{array}$$
Note that, by combining this system of differential equations into a single equation to eliminate ∂ρ∕∂r and ∂T∕∂r in Eq. (20), we can see that these equations still contain the same critical point where ${v}^{2}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{c}_{\text{s}}^{2}+{v}_{\text{turb}}^{2}$, although now with c_{s} a function of r.
In order to solve these differential equations, we begin by assuming initial velocity and density profiles $$\begin{array}{ll}v(r)\hfill & \text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{\sqrt{{c}_{\text{s}}^{2}+{v}_{\text{turb}}^{2}}}{1{R}_{\text{p},\text{mod}}/r}\left(1\frac{{R}_{\ast}}{r}\right)\hfill \\ \rho (r)\hfill & \text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{\dot{M}}{4\text{\hspace{0.17em}}\pi \text{\hspace{0.17em}}v(r)\text{\hspace{0.17em}}{r}^{2}}.\hfill \end{array}$$
As boundary conditions, we fix R_{*} and impose that T(R_{*}) = T_{eff}. These are then sufficient conditions to begin the following iteration procedure.
 1.
A radial grid is built up inward (toward R_{*}) and outward (away from R_{*}) from the critical point, R_{p,mod}, as determined in the prior iteration (or the initial conditions for the first iteration). Radial points are spaced by H(r)∕3 to properly resolve the variations in all hydrodynamic quantities.
 2.
The prior structures in velocity, density, and temperature are interpolated onto the new radial grid.
 3.
Gradients of all necessary quantities for Eqs. (19)–(21) are computed at each grid point. At the critical point, r = R_{p,mod}, where the velocity gradient is illdefined, l’Hopital’s rule is used.
 4.
Using these gradients, Eqs. (19)–(21) are numerically solved at each point on the radial grid inward and outward from the critical point, R_{p,mod}.
 5.
The resulting velocity, density, and temperature as functions of radius are used to begin the process again from Step 1.
This loop is repeated until velocity, density, and temperature at each point agree from one iteration to the next to better than 1%. At this point, the total optical depth of the wind is computed at the current mass loss rate. As the stellar radius is defined to be at τ(R_{*}) = 2∕3, the stellar massloss rate is recomputed according to $${\dot{M}}_{\text{new}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{\dot{M}}_{\text{old}}\frac{2/3}{{\tau}_{\text{old}}({R}_{\ast})}.$$(24)
As optical depth is proportional to density, which is in turn proportional to Ṁ, a constant change in massloss rate and thus density at every point has the effect of forcing τ(R_{*}) = 2∕3, assuming v(r) is unchanged. As changing ρ(r) will certainly change v(r), however, the final velocity profile is combined with the updated mass loss as new initial conditions to restart the loop described above. This outer iteration to update massloss rate is repeated until the update in Ṁ is also less than 1%, at which the final mass loss rate is returned as Ṁ_{num}. Note that for the remainder of the paper we include the subscript “num” to differentiate this numerically determined, nonisothermal massloss rate from the isothermal, analytic one derived in Sect. 2, which we henceforth call Ṁ_{an}.
Here itis again interesting to examine our “typical” model with M_{*} = 10 M_{⊙}, R_{*} = 500 R_{⊙}, T_{eff} = 3500 K, v_{turb} = 15 km s^{−1}, and κ = 0.01 cm^{2} g^{−1}. Using the method described above, we find that it has a nonisothermal massloss rate of Ṁ_{num} = 4.4 × 10^{−7}M_{⊙} yr^{−1}, approximately half the analytic rate. We can also take this opportunity to plot the radial profiles of density, velocity, and temperature as shown in Fig. 5. In all panels we use a vertical, red, dashed line to denote the numerically determined R_{p,mod}. As was the case in Sect. 2.1, optical depth is predominantly accumulated below R_{p,mod} such that the massloss rate in each simulation is primarily determined by integration over only these near star regions. However, like the isothermal, thermal pressuredriven Parker wind, the total driving energy available to the wind diverges as the outer radius goes to infinity because v_{turb} is constant, such that v_{∞} is not uniquely determined by such a model. While the terminal velocity of the wind will depend on the details of the radial v_{turb} profile, radial temperature profile, and any additional forces acting on the outer wind (for instance dust opacity), it is promising that these simulations already generate velocities at 30 R_{*} comparable to typically inferred values of v_{∞} for red supergiants.
While we here present the massloss rate from a simplified gas temperature structure arising from radiative equilibrium, in reality red supergiant winds may have regions that are heated by shocks and dissipation of turbulence (see, e.g. Schrijver & Zwaan 2000). In general, such shock heating and dissipation is not inferred to generate wind temperatures throughout the full wind acceleration region significantly in excess of the stellar effective temperature, nor should the wind be expected in general to cool below radiative equilibrium. Therefore, the simplified cases we present here of an isothermal wind at the stellar effective temperature and a wind with temperature set according to the spherically modified gray atmosphere can be roughly taken as thermal structures leading to, respectively, upper and lower limits to the massloss rate. For individual objects where more is known about the temperature profile, however, the method presented in this section can be directly applied by replacing the temperature information in Eqs. (17) and (21). As a sample case, let us take the chromospheric temperature profile inferred by O’Gorman et al. (2020) for Antares. In these observations, gas temperature is inferred to decrease from the stellar effective temperature to a minimum around 1.35R_{*} before coming back up to a maximum approximately equal to the stellar effective temperature around 2.5 R_{*}, and then finally falling off outward from there (see O’Gorman et al. 2020, their Fig. 3). Inserting this temperature profile into our model with all other parameters still as defined in our “typical” model above returns a massloss rate 7.7 × 10^{−7}M_{⊙} yr^{−1}. As anticipated, this falls between the limiting massloss rates from an isothermal model at the stellar effective temperature (8.5 × 10^{−7}M_{⊙} yr^{−1}) and a model using the Lucy (1971) spherically modified gray atmosphere (4.4 × 10^{−7}M_{⊙} yr^{−1}).
Fig. 5 Radial dependence of density (left), velocity (center), and temperature (right) for the nonisothermal case with M_{*} = 10 M_{⊙}, R_{*} = 500 R_{⊙}, T_{eff} = 3500 K, v_{turb} = 15 km s^{−1}, and κ = 0.01 cm^{2} g^{−1}. The red dashed line denotes the modified Parker radius. 
Fig. 6 Ratio of Ṁ calculated numerically for a nonisothermal situation to the analytic Ṁ calculated for an isothermal wind plotted as a function of v_{turb} and colored by M_{*}, varying R_{*}, M_{*}, and v_{turn}. The effectivetemperature and Rosseland opacity are held fixed at T_{eff} = 3500 K and κ_{R} = 0.01 cm^{2} g^{−1}. 
3.2 Nonisothermal correction factor
To investigate the role of the radial decline in temperature on the massloss rate, we run a grid of steadystatemodels varying stellar mass, stellar radius, and turbulent velocity and then plot the ratio Ṁ_{num}∕Ṁ_{an}. For this test we fix T_{eff} = 3500 K, motivated by the only modest dependence of Ṁ on this property (see Fig. 2), and continue to use κ_{R} = 0.01 cm^{2} g^{−1}. As shown by Fig. 6, this correction factor is itself a function of the model parameters that we can fit. We choose two methods of doing so, one where we fit the product of a constant normalization and powerlaws over stellar mass, stellar radius, and turbulent velocity individually (hereafter Method 1), and one where we apply our knowledge that the dominant competition in getting material off the star is between gravity and turbulent pressure to replace the powerlaw fits over each parameter independently with a single powerlaw fit over v_{turb}∕v_{esc}(M_{*}, R_{*}) (hereafter Method 2).
As figures of merit to compare the two fits, we use χ^{2}, the mean and standard deviation of the difference between the fit formula and the actual ratio Ṁ_{num}∕Ṁ_{an}, and the maximum error between the fit and Ṁ_{num}∕Ṁ_{an}. The comparison of the best models from each method is summarized in Table 1. Given that Method 1 uses four free parameters (a mean and three powerlaw indices) versus the two free parameters used by Method 2 (a mean and one powerlaw index), and that the figures of merit we have selected to compare between the two methods are comparable, we propose that the twoparameter fit from Method 2 should be preferable moving forward. Therefore, the nonisothermal correction factor is $$\left(\frac{{\dot{M}}_{\text{num}}}{{\dot{M}}_{\text{an}}}\right)\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{\left(\frac{{v}_{\text{turb}}/(17\text{\hspace{0.17em}}\text{km}\text{\hspace{0.17em}}{\text{s}}^{1})}{{v}_{\text{esc}}({M}_{\ast},{R}_{\ast})/(60\text{\hspace{0.17em}}\text{km}\text{\hspace{0.17em}}{\text{s}}^{1})}\right)}^{1.30},$$(25)
where we have omitted the leading constant as the best fit for this constant was 10^{0.00} = 1.0.
Note from Table 1 that this fit is quite good, with the correction factor recovered to about a factor of two in the worstcase for the two parameter fit. Therefore we present the combination of Eq. (25) with Eqs. (5), (13), and (14) as a predictive, theoretical massloss rate scaling for red supergiants.
Comparison of the fit methods.
4 Comparison of the theoretical massloss rate scaling with empirical rates and observed massloss
Now that we have this theoretical, predictive massloss rate as a function of stellar parameters, a natural next step is to compare this with empirical fits to red supergiant mass loss and individual observations. Given that no systematic theoretical predictions for RSG mass loss exist at the present, such empirical fits constitute the state of the art regarding RSG mass loss in applications such as stellar evolution. We review the empirical rates implemented in the stellar evolution code MESA, as well as the recently published rates from Beasor et al. (2020) in Sect. 4.1 in order to build up a representative sample. In Sect. 4.2 we then compare these empirical rates to our theoretical ones. Here we highlight that the empirical rates in general do not depend on all three fundamental stellar parameters (i.e. M_{*}, R_{*}, and L_{*}) while the theoretical massloss rates presented here does. Therefore Sect. 4.2 also includes a discussion of the scalings of stellar luminosity, mass, and radius with one another that we use to perform the comparison of our rates with the empirical ones. Finally, in Sect. 4.3 we compare our predicted rates with observationally inferred massloss rates for a selection of red supergiants.
4.1 Background on empirical rates used so far
The current version of the stellar evolution code MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019) provides builtin access to three empirical massloss rates for red supergiants. These come from de Jager et al. (1988) $$\dot{M}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{10}^{8.158}{\left(\frac{{L}_{\ast}}{{L}_{\odot}}\right)}^{1.769}{\left({T}_{\text{eff}}\right)}^{1.676},$$(26)
Nieuwenhuijzen & de Jager (1990) $$\dot{M}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{10}^{14.02}{\left(\frac{{L}_{\ast}}{{L}_{\odot}}\right)}^{1.24}{\left(\frac{{M}_{\ast}}{{M}_{\odot}}\right)}^{0.16}{\left(\frac{{R}_{\ast}}{{R}_{\odot}}\right)}^{0.81}$$(27)
and van Loon et al. (2005) $$\dot{M}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{10}^{5.65}{\left(\frac{{L}_{\ast}}{{10}^{4}{L}_{\odot}}\right)}^{1.05}{\left(\frac{{T}_{\text{eff}}}{3500}\right)}^{6.3}.$$(28)
To these we add the recently published empirical rate from Beasor et al. (2020) $$\dot{M}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{10}^{26.40.23{M}_{\ast}/{M}_{\odot}}{\left(\frac{{L}_{\ast}}{{L}_{\odot}}\right)}^{4.8},$$(29)
as the authors identify that this rate is notably lower than those previously published, a fact that we use to provide an indication to the range of possible observationally inferred rates. While these four are far from the only empirical massloss prescriptions that are used, they provide a representative sample. For a review of a variety of other empirical rates, see, for example, Mauron & Josselin (2011).
Fig. 7 Luminosity as a function of mass from red supergiants in Paxton et al. (2011) and Ekström et al. (2012). The overplotted curves assume ${L}_{\ast}/{L}_{\odot}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}f{({M}_{\ast}/{M}_{\odot})}^{3}$ with f = 12.5 for the MESA data and f = 18.5 for the Geneva data. 
4.2 Comparison of the theory model to empirical rates
As mentioned above, it is common for empirical rates to be provided as functions that do not include a dependency on all three fundamental stellar parameters (i.e. M_{*}, R_{*}, and L_{*}, or equivalently M_{*}, L_{*}, and T_{eff}). For the theoretical rates, however, all three parameters are needed. To obtain these, we consider stellar evolution tracks for nonrotating stars from both MESA (Paxton et al. 2011) and GENEC (Ekström et al. 2012). Note that GENEC uses a somewhat different prescription for the massloss rate than what is discussed above for MESA. For M_{*} ≤ 12 M_{⊙}GENEC adopts the Reimers (1975, 1977) rates, whereas for M_{*}≥ 15 M_{⊙}GENEC uses de Jager et al. (1988) for log_{10}(T_{eff}) > 3.7 and a linear fit to data from Sylvester et al. (1998) and van Loon et al. (1999) otherwise. Despite these differences in massloss prescription, comparison of Paxton et al. (2011) and Ekström et al. (2012) shows that in both cases the red supergiant branch roughly follows ${L}_{\ast}/{L}_{\odot}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}f{({M}_{\ast}/{M}_{\odot})}^{3}$ as shown by Fig. 7 with f = 12.5 for MESA and f = 18.5 for Geneva. Therefore, as a compromise we chose f = 15.5 for our massluminosity relation to compare with the empirical rates. As a simplifying assumption, we continue to hold T_{eff} fixed at 3500 K as was done in Sect. 3, so that ${L}_{\ast}\propto {R}_{\ast}^{2}$.
Using these massradiusluminosity relationships, we can now compare the empirical and theoretical rates as shown by Fig. 8. As the empirical rates differ from one another by several orders of magnitude at all luminosities, we select four values of v_{turb} from 15 to 21 km s^{−1} for this comparison. Despite this scatter between individual empirical rates and in turn between the theoretical rates using different turbulent velocities, the theoretical rates are broadly consistent with the empirical ones for reasonable values of v_{turb} (see, e.g. Josselin & Plez 2007; Ohnaka et al. 2017). Further examining this comparison, one notes that the empirical rates have shallower slopes than the theoretical rates. This may suggest that there is a trend of decreasing v_{turb} for more massive, more luminous, and larger radius red supergiants. Alternatively, as is suggested by the mismatch in both value and slope of the empirical rates with one another, this shallower slope may simply reflect some missing physics when inferring and calibrating these empirical rates or potential biases in the observed sample.
We note that the massluminosityradius relations that we use here may not be appropriate to Galactic stars starting out their lives rapidly rotating, as these tend to evolve to higher luminosities for the same initial mass (see, e.g. Heger & Langer 2000; Meynet & Maeder 2000; Brott et al. 2011). Given, however, the very similar effective temperatures of these stars on the red supergiant branch this is consistent with these stars evolving to have larger radii, lower surface gravities, and higher massloss rates if all other physics remains the same. This further suggests that two stars with different masses can reach the same point on the HR diagram with different massloss rates which would confound any massloss rate predictions omitting a dependence on M_{*}. In fact the general treatment that massloss rates are enhanced by rotation would go in the direction of only further exacerbating this discrepancy (see, e.g. Friend & Abbott 1986; Langer 1997; Heger & Langer 1998; Maeder & Meynet 2000). However, as further discussed in Sect. 5, these predictions of increased massloss rate with increased rotation rate are based on fits to the theoretical predictions of Friend & Abbott (1986). This work addresses only the wind massloss rate of hot stars with radiation driven outflows, and has been questioned by subsequent work (e.g. Owocki et al. 1996; Petrenz & Puls 2000) such that the question of the role of rotation in the massloss of red supergiants remains open.
Fig. 8 An overplot of three of the theoretical massloss rate scalings used in MESA with the predicted massloss rate from the theoretical scaling. Four values of v_{turb} are used as noted in the plot label. 
4.3 How does the theory model compare to what is observationally inferred?
While the preceding comparison of the empirical rates with theoretical ones helps to generally motivate the quality of the theoretical massloss rate predictions, we can go one step further and compare theory and observation for a sample of observed stars. For this we select the observations compiled by Josselin & Plez (2007) with the addition of Antares using data from Ohnaka et al. (2017). Note that this limited sample is based on these authors compiling velocity dispersion measurements that we can compare with the v_{turb} required by our model. Table 2 recapitulates the stellar masses, effective temperatures, and radii from Josselin & Plez (2007) and Ohnaka et al. (2017). As the model we present predicts gas massloss rates while Josselin & Plez (2007) report dust massloss rates, we apply a constant gas to dust mass ratio Ṁ_{gas} = 250 Ṁ_{dust}. This ratio is itself not particularly well constrained by the existing literature as shown by a sample comparison of the gas massloss rates of De Beck et al. (2010) with the dust massloss rates of Harwit et al. (2001) and Verhoelst et al. (2009) for VY CMa, μ Cep, and α Ori, which returns ratios between ~70 and ~700. Finally, we report the velocity dispersion measurement that Josselin & Plez (2007) and Ohnaka et al. (2017) associate with turbulent velocity as v_{turb, Obs}, and the turbulent velocity that would be required for the theoretical model to reproduce Ṁ_{gas} as v_{turb, Theory}. For the data from Josselin & Plez (2007), we round their reported values to the nearest km s^{−1} as the errors on the original data may beas large as ±2 km s^{−1}. We report our requiredv_{turb} with the same accuracy.
Comparing the mean value of the observed turbulent velocity ${\overline{v}}_{\text{turb},\text{\hspace{0.05em}}\text{Obs}}$ to the mean theoretical turbulent velocity ${\overline{v}}_{\text{turb},\text{\hspace{0.05em}}\text{Theory}}$, the two are consistent with each other over the population to within one standard deviation, as ${\overline{v}}_{\text{turb},\text{\hspace{0.05em}}\text{Obs}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}20.2\text{\hspace{0.17em}}\pm \text{\hspace{0.17em}}3.4$ km s^{−1} and ${\overline{v}}_{\text{turb},\text{\hspace{0.05em}}\text{Theory}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}18.2\text{\hspace{0.17em}}\pm \text{\hspace{0.17em}}3.5$ km s^{−1}. Repeating this exercise with the extremal gas to dust mass ratios above also retrieves ${\overline{v}}_{\text{turb}}$ values consistent with the ${\overline{v}}_{\text{turb},\text{\hspace{0.05em}}\text{Obs}}$ and ${\overline{v}}_{\text{turb},\text{\hspace{0.05em}}\text{Theory}}$ values reported here to within one standard deviation, such that this choice has not significantly impacted our results. These ranges are also generally consistent with the turbulent velocities required to recover the empirical rates in the prior subsection. Moreover, the overall ranges are consistent with v_{turb, Obs} ∈ [12, 24] km s^{−1} and v_{turb, Theory} ∈ [13, 26] km s^{−1}. Finally, looking at the distributions of v_{turb} for both theory and observation, plotted in the left panel of Fig. 9, as well as comparing v_{turb, Obs} to v_{turb, Theory} for individualobjects as plotted in the right panel of Fig. 9, one can see that the majority of turbulent velocities required by theory (10 out of 13) are lower than what is observationally inferred. This argues that turbulent pressure is an ample driving source for the wind of red supergiants, strongly reinforcing the quality of this model in explaining the general behavior of these objects’ mass loss.
As turbulent velocity is not observationally inferred for many red supergiants, v_{turb} = 18.2 km s^{−1} as the mean v_{turb, Theory} found above may thus bea reasonable value to use for application of our model in future works. Similarly, a reasonable range of variations on this would be 14.8 km s^{−1} < v_{turb} < 21.6 km s^{−1} as suggested by plus or minus one standard deviation about this mean. While this will not provide a perfect match to the massloss rates of individual objects (see for example the right panel of Fig. 9) it provides a representative prediction for the average massloss behavior of red supergiants such that it would be appropriate to use in stellar evolution codes. As such, and as discussed further in the following section, future work should be dedicated to improving this estimation, establishing whether v_{turb} varies in time for individual stars, and testing its dependence on key stellar parameters.
Observed parameters from Josselin & Plez (2007) with the addition of Antares as observed by Ohnaka et al. (2013, 2017), including v_{turb} as required by theory to recover the massloss rate keeping all other parameters fixed.
5 Summary and future directions
We have here derived a theoretical, analytic (dustfree) model for the mass lossrates of red supergiants. Building upon the works by Gustafsson & Plez (1992) and Josselin & Plez (2007), we examine the role of turbulent pressure in levitating material in the stellar atmosphere up to a modified equivalent of the standard Parker wind critical point. This is achieved by employing a standard decomposition of pressure into thermal and turbulent pressure components as, for instance, employed in standard 1D spherically symmetric model atmospheres (e.g. Gustafsson et al. 2008). By applying the resulting expression for pressure in the steadystate equations for conservation of mass and momentum, we examine both hydrostatic and low velocity levitation for an isothermal model (Sect. 2). We also numerically solve the full 1D steadystate equationofmotion for a nonisothermal model employing the Lucy (1971) temperature structure as a stand in for an energy equation. By combining each of these cases of steadystate levitation with the constraint that the total accumulated optical depth of circumstellar material from infinity down to the stellar photosphere at R_{*} = r(τ = 2∕3), we provide unique massloss rates.
As shown in Sect. 2 this method provides a fully analytic massloss rate for an isothermal wind. Using a grid of numerical models, we further provide a fit to the correction factor between the isothermal and nonisothermal massloss rate for the same input parameters (Sect. 3.2). Combining this correction factor with the isothermal massloss rates then presents a fully analytic massloss rate prescription for red supergiants, given by the combination of Eqs. (5), (13), (14), and (25). Comparisons of this massloss rate with both observations and empirical massloss rate prescriptions show that such a turbulent pressure driven massloss rate is able to recover realistic massloss rates over the range of parameters appropriate to red supergiants. These comparisons also inform our suggested choice of v_{turb} = 18.2 ± 3.4 km s^{−1} as an appropriate representative value of v_{turb} in massloss rate calculations where this parameter is not well constrained.
To continue the discussion, we enumerate some items of future work required to place this analytic, theoretical model of red supergiant massloss on even firmer footing. From an observational side, as an immediate first item estimations of turbulent velocities for a wider body of red supergiants are needed. As it stands, the limited sample that we have examined here provides highly encouraging results, and additional objects would allow the model to be calibrated more carefully. Additionally, such observations would be of key importance to understand how this current free parameter of the model scales with stellar parameters. Finally, such future observations could also help constrain the depth dependence of v_{turb}. Prior work has suggested that turbulent velocity may actually increase with distance from the star over at least part of the wind launching region (de Koter et al. 1988; Josselin & Plez 2007), and better understanding this profile would be highly useful in further developing the model presented here.
To further develop the model itself, main efforts may focus on relaxing key approximations made in this paper. One of the first of these would be to use a radial profile of v_{turb}, perhaps as inferred by observations above, rather than the spatial average employed here. Such a radial profile might arise from dissipation of turbulence by shocks. Inclusion of a radial profile in v_{turb} would allow,for instance, for a prediction of the terminal velocities of red supergiant winds instead of only massloss rates, as the energy available to drive the wind would no longer diverge with increasing outer radius of the model. Further generalization of v_{turb} could also include variations over the surface of the star and/or in time. Both of these would arise naturally if turbulence is related to surface convection in the star, and their inclusion would help to explain the observed clumpy structure of the wind.
A further investigation of the overall structure of red supergiant winds may also reintroduce dust and/or molecules and their associated opacity to examine the impact this would have for the proposed model. As discussed in the introduction, previous models have assumed that radiation pressure is the dominant driving force in carrying the mass loss of red supergiants, with the main question being the origin of the “missing force” necessary to get material from the stellar photosphere up to the region where molecules and dust can form, the latter occurring potentially several stellar radii away from the photosphere. Seemingly circumventing such a model, we show that turbulent pressure extended red supergiant atmospheres are already able to generate the appropriate massloss rates in the absence of radiation pressure on dust or molecules. However, it is important to recognize that this does not necessarily mean that enhanced opacity from dust or molecules has a negligible contribution. In fact, these opacities almost certainly play a crucial role in setting the wind terminal speed as turbulent pressure dies off away from the stellar surface (see the discussion in Sect. 3.1). It is also clear, as discussed in Sect. 2.2, that the overall scale of opacity and its radial profile may contribute to the optical depth and scale height of the wind and thus to the mass loss rate retrieved by this model. Further, it is even possible that the radial profile of turbulent pressure is such that turbulent pressure plays the dominant role in levitating material off the stellar surface, but that dust or molecular opacity carries the wind through the critical point as previously theorized.
Another important line of investigation regards the origin of the turbulent velocities and extended atmospheres this model leverages. As discussed by ArroyoTorres et al. (2015), previous models have focused on convection through 3D, socalled “starinabox” radiation hydrodynamics simulations such as those computed with CO^{5}BOLD (Freytag et al. 2012), and on pulsations through 1D simulations akin to those used for AGB stars and especially Mira variabls (e.g. CODEX, Ireland et al. 2008, 2011). While convection and pulsations in the stellar atmosphere both provide tempting options, convection simulations generally fail to generate the turbulent velocities or atmospheric extensions inferred from observations, and pulsations are required to be of unrealistically high amplitude compared to observed lightcurves for red supergiants (see, e.g. Josselin & Plez 2007; ArroyoTorres et al. 2015).
Similarly, one can examine the role that rotation might have in generating similar levels of extension to the stellar atmosphere. This can be done by applying the centrifugal term (Ω^{2} r) either along side or instead of v_{turb} in the preceding analysis. Taking the limiting case where turbulent pressure is omitted, we find that matching the atmospheric extension of a 18.2 km s^{−1} turbulent velocity even at the most optimal position on the stellar equator (including the effects of stellar oblateness) requires a rotation velocity almost 94% of the orbital velocity, implausibly fast for a red supergiant. As we reinforce through this paper that atmospheric extension is crucial for the massloss process, further examination of the possible missing ingredients for such models as well as potential constructive interactions between them is now highly timely.
To wrap up the discussion, we address some of the potential implications of applying our massloss rates, especially for the case of stellar evolution modeling. As discussed in Sect. 4.2, the current stateoftheart massloss rates for stellar evolution modeling are all empirical. Given that none of the empirical rates match the theoretical one we present here for all stars, these new theoretical rates should impact the mass, luminosity, and radius distributions of supergiant supernova progenitors. Additionally, the steeper trend of increasing massloss rate with increasing luminosity would naturally turn more massive stars evolving toward the red supergiant branch back to the blue due to envelope stripping, thereby predicting a decreased upper luminosity limit for red supergiants. Finally, these rates provide intriguing implications for stellar evolution in lower metallicity environments. Many standard stellar evolution implementations of mass loss in the red supergiant branch impose a downturn with metallicity partially inspired by the Ṁ ∝ Z^{n} dependence of linedriven hot star winds (Vink et al. 2001; Mokiem et al. 2007; Björklund et al. 2020). Observations analyzed by Mauron & Josselin (2011) do seem to suggest such a scaling of $\dot{M}\propto {\left(Z/{Z}_{\odot}\right)}^{0.7}$ from the Galaxy to the red supergiants in the Small Magellanic Cloud. However, they concluded that the scaling to the Large Magellanic Cloud was not well constrained by the samples they considered, emphasizing that more investigation is needed. Meanwhile, the model we present here should not be strongly dependent on metallicity. While the optical depth and scale height of the wind clearly do depend on dust and molecular content, which in turn depend on metallicity, this is a weaker secondary scaling when compared with the turbulent pressure itself. Moreover, decreasing opacity can lead to increased massloss rates in this model depending where in parameter space a star sits (see Sect. 2.2). So long as the turbulent pressure mechanism itself is not strongly metallicity dependant, which turbulence triggered by hydrogen recombination is not, this could then imply much stronger mass loss on the red supergiant branch at low metallicity than previously assumed, and as such significantly different stellar evolution tracks.
Fig. 9 Comparisons of v_{turb} from observations (v_{turb, Obs}) and frommatching corresponding observationally inferred massloss rates to the theoretical predictions presented here (see text). Left panel: histogram of the distribution while right panel: plots v_{turb, Obs} and v_{turb, Theory} for each object numbered according to Table 2. 
Acknowledgements
We would like to thank the anonymous referee for their comments which enhanced the discussion in this paper. We acknowledge support from the KU Leuven C1 grant MAESTRO C16/17/007. L.D. also acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreements No. 646758: AEROSOL with PI L. Decin).
Appendix A Computation of g_{rad} and κ from accelerating molecules
As discussed in Sect. 2, the total opacity scale of the winds we treat here can be strongly impacted by the Doppler shifting of spectral lines, leading to opacities exceeding the static Rosseland mean opacity. To show this, we begin with a generalized form of radiative acceleration in the $\widehat{n}$ direction, $${g}_{\text{rad}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{1}{c}{\displaystyle \oint {\displaystyle {\int}_{0}^{\infty}{\kappa}_{\nu}}}\text{\hspace{0.17em}}{I}_{\nu}\text{\hspace{0.17em}}\text{d}\nu \text{\hspace{0.17em}}\widehat{n}\text{\hspace{0.17em}}\text{d}\Omega ,$$(A.1)
with speed of light c, extinction coefficient κ_{ν} [cm^{2} g^{−1}], and intensity per unit frequency I_{ν}. If one assumes that the extinction comes from a single, isolated spectral line, then a natural choice is to split κ_{ν} into a normalized, frequencydependent shape of the spectral line, or profile function ϕ_{ν}, and a lineintegrated total extinction coefficient $${\kappa}_{L}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{{\sigma}_{\text{cl}}\text{\hspace{0.17em}}{n}_{\text{l}}\text{\hspace{0.17em}}{f}_{\text{lu}}}{\rho}\left(1\mathrm{exp}\left(\frac{h\text{\hspace{0.17em}}{\nu}_{\text{o}}}{{k}_{\text{B}}\text{\hspace{0.17em}}T}\right)\right).$$(A.2)
In this expression for κ_{L}, σ_{cl} is the classical oscillator cross section, n_{l} is the number density of particles in the lower level of the transition, f_{lu} is the quantum mechanical oscillator strength of the transition, and ρ is mass density. The final exponential term accounts for stimulated emission under the assumption of local thermodynamic equilibrium (LTE). Such stimulated emission can play an important role for the considered infrared spectral lines transitions as the energy of the transition, given by Planck’s constant h times rest frequency ν_{o}, can be comparable to the thermal energy of the gas, given by Boltzmann’s constant k_{B} times gas temperature T. Note that here κ_{L} has picked up units of [cm^{2} g^{−1} Hz] as the profile function is defined per unit frequency. Under these substitutions, the equation we now need to solve becomes $${g}_{\text{line}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{{\kappa}_{L}}{c}{\displaystyle \oint {\displaystyle {\int}_{0}^{\infty}{\varphi}_{\nu}}}\text{\hspace{0.17em}}{I}_{\nu}\text{\hspace{0.17em}}\text{d}\nu \text{\hspace{0.17em}}\widehat{n}\text{\hspace{0.17em}}\text{d}\Omega .$$(A.3)
At this point, it is important to recognize that the intensity fed into this expression can itself be a strong function of frequency even in the presence of an optically thin continuum, as is the case for red supergiants, due to the attenuation of intensity by the spectral line we consider itself. Therefore, we have to obtain intensity from the formal solution of the radiative transfer equation $${I}_{\nu}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{I}_{\nu}^{\text{o}}{\text{e}}^{{\tau}_{\nu}}+{\displaystyle \int {S}_{\nu}}({t}_{\nu})\text{\hspace{0.17em}}{\text{e}}^{\text{}{\tau}_{\nu}{t}_{\nu}\text{}}\text{\hspace{0.17em}}\text{d}{t}_{\nu}.$$(A.4)
While in general this becomes quite complex as the source function S_{ν} depends on I_{ν}, we can simplify by leveraging symmetry arguments. In the Sobolev limit considered here (see further below), the second term here becomes foreaft symmetric such that it cancels out in the integral $\oint \widehat{n}}\text{d}\Omega $. As regards the optical depth τ_{ν}, we take the approximation that this arises only from extinction from the spectral line itself such that $${\tau}_{\nu}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{\displaystyle \int {\kappa}_{L}}\text{\hspace{0.17em}}{\varphi}_{\nu}\text{\hspace{0.17em}}\rho \text{\hspace{0.17em}}\text{d}l,$$(A.5)
where the mixed frequency and spatial notations can be rectified by using the Doppler formula $$\frac{\nu {\nu}_{\text{o}}}{{\nu}_{\text{o}}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{{v}_{l}}{c},$$(A.6)
such that $${\tau}_{\nu}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{\displaystyle \int {\kappa}_{L}}\text{\hspace{0.17em}}{\varphi}_{\nu}\text{\hspace{0.17em}}\rho \left(\frac{\partial l}{\partial {v}_{\text{l}}}\right)\left(\frac{\partial {v}_{\text{l}}}{\partial \nu}\right)\text{d}\nu \text{\hspace{0.17em}}=\text{\hspace{0.17em}}{\displaystyle \int \frac{{\kappa}_{L}\text{\hspace{0.17em}}\rho \text{\hspace{0.17em}}c}{{\nu}_{\text{o}}\text{\hspace{0.17em}}\partial {v}_{l}/\partial l}}\text{\hspace{0.17em}}{\varphi}_{\nu}\text{\hspace{0.17em}}\text{d}\nu .$$(A.7)
Finally, we can take the Sobolev approximation (Sobolev 1960) which argues that all hydrodynamic variables are effectively constant over a line resonance region and therefore can be pulled through the frequency integral to give $${g}_{\text{line}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{{\kappa}_{L}}{c}{\displaystyle \oint {\displaystyle {\int}_{0}^{\infty}{\varphi}_{\nu}}}\text{\hspace{0.17em}}{I}_{\nu}^{\text{o}}\text{\hspace{0.17em}}{\text{e}}^{{\tau}_{S}{\displaystyle {\int}_{0}^{\infty}{\varphi}_{\nu}}\text{\hspace{0.05em}}\text{d}\nu}\text{d}\nu \text{\hspace{0.17em}}\widehat{n}\text{\hspace{0.17em}}\text{d}\Omega ,$$(A.8)
where we have introduced the Sobolev optical depth τ_{S} = κ_{L}ρc∕(ν_{o}∂v_{l}∕∂l). This final expression can now be analytically solved to give $${g}_{\text{line}}=\frac{{\kappa}_{L}}{c}{\displaystyle \oint {I}_{\nu}^{\text{o}}}\left(\frac{1{e}^{{\tau}_{S}}}{{\tau}_{S}}\right)\text{d}\Omega ,$$(A.9)
which further simplifies in the case of a spherically symmetric flow and radiation from a point star to $${g}_{\text{line},\text{r}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{\pi \text{\hspace{0.17em}}{\kappa}_{L}\text{\hspace{0.17em}}{I}_{\nu}^{\text{o}}}{c}\left(\frac{1{\text{e}}^{{\tau}_{\text{S},\text{r}}}}{{\tau}_{\text{S},\text{r}}}\right),$$(A.10)
with τ_{S,r} replacing the general line of sight velocity gradient ∂v_{l}∕∂l in τ_{S} with the radial velocity gradient ∂v_{r}∕∂r.
The Sobolev approximation which we have employed here is valid in wind outflow regions where the flow speed exceeds a few times the characteristic thermal speed, v_{th}, such that there is not significant absorption of intensity entering the resonance region by a spectral line in the hydrostatic stellar atmosphere, and such that the physical size of the resonance region is small compared to the scale over which the hydrodynamic quantities entering τ_{S} vary. In practice, this condition is nearly always met, however, as the thermal speeds of the molecules we examine arereduced compared to the bulk wind sound speed by the potentially quite substantial ratio of mean molecular mass of the bulk gas to mass of the molecule. Therefore, as this expression holds for each individual line over the bulk of the wind, to generalize this to a spectrum of lines one simply needs to sum g_{line}. Numerous tabulated line lists are available to facilitate this process, each including ν_{o}, either f_{lu} or the Einstein coefficients necessary to compute it, and the energy and degeneracy of the levels involved in each of the transitions, which can be used to compute n_{l}, for instance from the Boltzmann distribution in Local Thermodynamic Equilibrium. Taking this fully summed acceleration, it is then straightforward to back out a fluxweighted mean opacity κ of the line acceleration by appealing to Eq. (3) such that $$\kappa \text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{c}{{F}_{\text{r}}}{\displaystyle \sum {g}_{\text{line},\text{r}}},$$(A.11)
where the summation is over all lines. As more force is available from an accelerating spectral line than a static one, this will always return fluxweighted mean opacity that is larger than the Rosseland mean opacity, thereby substantiating our choice to test the effects of varying κ in Sect. 2.2.
References
 ArroyoTorres, B., Wittkowski, M., Chiavassa, A., et al. 2015, A&A, 575, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Barber, R. J., Tennyson, J., Harris, G. J., & Tolchenov, R. N. 2006, MNRAS, 368, 1087 [NASA ADS] [CrossRef] [Google Scholar]
 Beasor, E. R., Davies, B., Smith, N., et al. 2020, MNRAS, 492, 5994 [CrossRef] [Google Scholar]
 Björklund, R., Sundqvist, J. O., Puls, J., & Najarro, F. 2020, A&A, accepted [arXiv:2008.06066] [Google Scholar]
 Bladh, S., & Höfner, S. 2012, A&A, 546, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bladh, S., Höfner, S., Nowotny, W., Aringer, B., & Eriksson, K. 2013, A&A, 553, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Braun, K., Baade, R., Reimers, D., & Hagen, H. J. 2012, A&A, 546, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 De Beck, E., Decin, L., de Koter, A., et al. 2010, A&A, 523, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Decin, L., Hony, S., de Koter, A., et al. 2006, A&A, 456, 549 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 de Jager, C., Nieuwenhuijzen, H., & van der Hucht, K. A. 1988, A&AS, 72, 259 [NASA ADS] [Google Scholar]
 de Koter, A., de Jager, C., & Nieuwenhuijzen, H. 1988, A&A, 200, 146 [Google Scholar]
 Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Freytag, B., Steffen, M., Ludwig, H. G., et al. 2012, J. Comput. Phys., 231, 919 [Google Scholar]
 Friend, D. B., & Abbott, D. C. 1986, ApJ, 311, 701 [NASA ADS] [CrossRef] [Google Scholar]
 Goorvitch, D. 1994, ApJS, 95, 535 [NASA ADS] [CrossRef] [Google Scholar]
 Gustafsson, B., & Plez, B. 1992, Instabilities in Evolved Super and Hypergiants, eds. C. de Jager, & H. Nieuwenhuijzen (Royal Netherlands Academy), 86 [Google Scholar]
 Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Harwit, M., Malfait, K., Decin, L., et al. 2001, ApJ, 557, 844 [NASA ADS] [CrossRef] [Google Scholar]
 Hearn, A. G. 1975, A&A, 40, 355 [NASA ADS] [Google Scholar]
 Heger, A., & Langer, N. 1998, A&A, 334, 210 [NASA ADS] [Google Scholar]
 Heger, A., & Langer, N. 2000, ApJ, 544, 1016 [NASA ADS] [CrossRef] [Google Scholar]
 Höfner, S., & Olofsson, H. 2018, A&ARv, 26, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
 Ireland, M. J., Scholz, M., & Wood, P. R. 2008, MNRAS, 391, 1994 [NASA ADS] [CrossRef] [Google Scholar]
 Ireland, M. J., Scholz, M., & Wood, P. R. 2011, MNRAS, 418, 114 [NASA ADS] [CrossRef] [Google Scholar]
 Josselin, E., & Plez, B. 2007, A&A, 469, 671 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Langer, N. 1997, ASP Conf. Ser., 120, 83 [Google Scholar]
 Levesque, E. M. 2017, Astrophysics of Red Supergiants (IOP Publishing) [CrossRef] [Google Scholar]
 Lucy, L. B. 1971, ApJ, 163, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Maeder, A., & Meynet, G. 2000, A&A, 361, 159 [NASA ADS] [Google Scholar]
 Mauron, N., & Josselin, E. 2011, A&A, 526, A156 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 McKemmish, L. K., Masseron, T., Hoeijmakers, H. J., et al. 2019, MNRAS, 488, 2836 [Google Scholar]
 Meynet, G., & Maeder, A. 2000, A&A, 361, 101 [NASA ADS] [Google Scholar]
 Mokiem, M. R., de Koter, A., Vink, J. S., et al. 2007, A&A, 473, 603 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nieuwenhuijzen, H., & de Jager, C. 1990, A&A, 231, 134 [NASA ADS] [Google Scholar]
 O’Gorman, E., Harper, G. M., Ohnaka, K., et al. 2020, A&A, 638, A65 [CrossRef] [EDP Sciences] [Google Scholar]
 Ohnaka, K., Hofmann, K. H., Schertl, D., et al. 2013, A&A, 555, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ohnaka, K., Weigelt, G., & Hofmann, K. H. 2017, Nature, 548, 310 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P., Cranmer, S. R., & Gayley, K. G. 1996, ApJ, 472, L115 [NASA ADS] [CrossRef] [Google Scholar]
 Parker, E. N. 1958, ApJ, 128, 664 [Google Scholar]
 Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
 Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
 Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34 [Google Scholar]
 Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10 [Google Scholar]
 Petrenz, P., & Puls, J. 2000, A&A, 358, 956 [NASA ADS] [Google Scholar]
 Reimers, D. 1975, Mem. Soc. R. Sci. Liege, 8, 369 [NASA ADS] [Google Scholar]
 Reimers, D. 1977, A&A, 61, 217 [NASA ADS] [Google Scholar]
 Rogerson, J. B., J., & Lamers, H. J. G. L. M. 1975, Nature, 256, 190 [Google Scholar]
 Schrijver, C. J., & Zwaan, C. 2000, Solar and Stellar Magnetic Activity (Cambridge: Cambridge University Press) [Google Scholar]
 Sobolev, V. V. 1960, Moving Envelopes of Stars (Cambridge: Harvard University Press) [Google Scholar]
 Sylvester, R. J., Skinner, C. J., & Barlow, M. J. 1998, MNRAS, 301, 1083 [NASA ADS] [CrossRef] [Google Scholar]
 van Loon, J. T., Groenewegen, M. A. T., de Koter, A., et al. 1999, A&A, 351, 559 [NASA ADS] [Google Scholar]
 van Loon, J. T., Cioni, M. R. L., Zijlstra, A. A., & Loup, C. 2005, A&A, 438, 273 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Verhoelst, T., van der Zypen, N., Hony, S., et al. 2009, A&A, 498, 127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wood, P. R., Bessell, M. S., & Fox, M. W. 1983, ApJ, 272, 99 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Observed parameters from Josselin & Plez (2007) with the addition of Antares as observed by Ohnaka et al. (2013, 2017), including v_{turb} as required by theory to recover the massloss rate keeping all other parameters fixed.
All Figures
Fig. 1 Optical depth calculated inward to r. Note that most of the optical depth is accumulated within a stellar radius of the star. 

In the text 
Fig. 2 Dependence of the analytic massloss rate on stellar mass (left), stellar radius (center), and temperature (right). 

In the text 
Fig. 3 Dependence of the analytic massloss rate on turbulent velocity. 

In the text 
Fig. 4 Dependence of the analytic massloss rate on Γ for M_{*} = 10 M_{⊙} and T_{eff} = 3500 K. Left panel: uses R_{*} = 500 R_{⊙} and v_{turb} = 10 km s^{−1}, middle panel: uses R_{*} = 500 R_{⊙} and v_{turb} = 15 km s^{−1}, and right panel: uses R_{*} = 1000 R_{⊙} and v_{turb} = 15 km s^{−1}. 

In the text 
Fig. 5 Radial dependence of density (left), velocity (center), and temperature (right) for the nonisothermal case with M_{*} = 10 M_{⊙}, R_{*} = 500 R_{⊙}, T_{eff} = 3500 K, v_{turb} = 15 km s^{−1}, and κ = 0.01 cm^{2} g^{−1}. The red dashed line denotes the modified Parker radius. 

In the text 
Fig. 6 Ratio of Ṁ calculated numerically for a nonisothermal situation to the analytic Ṁ calculated for an isothermal wind plotted as a function of v_{turb} and colored by M_{*}, varying R_{*}, M_{*}, and v_{turn}. The effectivetemperature and Rosseland opacity are held fixed at T_{eff} = 3500 K and κ_{R} = 0.01 cm^{2} g^{−1}. 

In the text 
Fig. 7 Luminosity as a function of mass from red supergiants in Paxton et al. (2011) and Ekström et al. (2012). The overplotted curves assume ${L}_{\ast}/{L}_{\odot}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}f{({M}_{\ast}/{M}_{\odot})}^{3}$ with f = 12.5 for the MESA data and f = 18.5 for the Geneva data. 

In the text 
Fig. 8 An overplot of three of the theoretical massloss rate scalings used in MESA with the predicted massloss rate from the theoretical scaling. Four values of v_{turb} are used as noted in the plot label. 

In the text 
Fig. 9 Comparisons of v_{turb} from observations (v_{turb, Obs}) and frommatching corresponding observationally inferred massloss rates to the theoretical predictions presented here (see text). Left panel: histogram of the distribution while right panel: plots v_{turb, Obs} and v_{turb, Theory} for each object numbered according to Table 2. 

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