Analytic, dust-independent mass-loss rates for red supergiant winds initiated by turbulent pressure
AAstronomy & Astrophysics manuscript no. rsg_modParker_final ©ESO 2021Monday 11 th January, 2021
Analytic, dust-independent mass-loss rates for red supergiantwinds initiated by turbulent pressure
N. D. Kee , J. O. Sundqvist , L. Decin , A. de Koter , , and H. Sana Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, B-3001 Leuven, Belgiume-mail: [email protected] Anton Pannekoek Institute, University of Amsterdam, Science Park 904, 1098XH Amsterdam, The Netherlands
ABSTRACT
Context.
Red supergiants are observed to undergo vigorous mass-loss. However, to date, no theoretical model hassucceeded in explaining the origins of these objects’ winds. This strongly limits our understanding of red supergiantevolution and Type II-P and II-L supernova progenitor properties.
Aims.
We examine the role that vigorous atmospheric turbulence may play in initiating and determining the mass-lossrates of red supergiant stars.
Methods.
We analytically and numerically solve the equations of conservation of mass and momentum, which we latercouple to an atmospheric temperature structure, to obtain theoretically motivated mass-loss rates. We then comparethese to state-of-the-art empirical mass-loss rate scaling formulae as well as observationally inferred mass-loss rates ofred supergiants.
Results.
We find that the pressure due to the characteristic turbulent velocities inferred for red supergiants is sufficientto explain the mass-loss 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 mass-loss rate prescription for redsupergiants. We conclude by highlighting some intriguing possible implications of these rates for future studies of stellarevolution, especially in light of the lack of a direct dependence on metallicity.
Key words.
Stars: mass-loss; Stars: winds, outflows; Stars: massive; Stars: supergiants; Turbulence
1. Introduction
The lack of a satisfactory theory explaining the strong, > − M (cid:12) yr − , mass loss for evolved massive stars onthe red supergiant branch has been a long-standing problemin our understanding of these objects (see Levesque 2017,for a recent discussion of this and other current puzzles instudying red supergiants). While this considerable gap inour knowledge has been patched over somewhat by empiri-cal rates and scaling formulae (see, e.g., Mauron & Josselin2011, for a review), the overall disagreement between anygiven two of these leaves those attempting to model redsupergiants directly, or attempting to model stellar evo-lution including or approaching this phase, with a ratheruntenable problem. Namely, one must somehow distinguishand select between rates that may differ from each otherby as much as four orders of magnitude in some parts ofthe parameter space without a general understanding of theunderlying process that determines this mass loss.In order to understand why theories generally struggleso much in this region, it is important to highlight somegeneral features of previous modeling attempts. In analogyto the lower-mass asymptotic giant branch (AGB) stars,it has generally been assumed that the primary drivingmechanism of red supergiant winds is radiation pressureon dust grains (see, e.g., Höfner & Olofsson 2018, for a re-view). For AGB stars, pulsations provide an atmosphericpiston levitating material off the stellar surface and up to a region where the temperature has dropped far enoughto allow dust to condense (e.g., Bladh & Höfner 2012;Bladh et al. 2013). However, as dust nucleation is also adensity-dependent process, it is essential not only that ma-terial reaches this region but that enough material is levi-tated to make the dust formation efficient. Due to both thelower pulsational amplitudes inferred for red supergiants(e.g., Wood et al. 1983) and their higher effective temper-atures relative to AGB stars, applying similar models tored supergiants thus requires gas to reach a greater heightin comparison to the stellar radius, while less material isexpected to be levitated in the first place. This implies astrongly decreased efficiency of dust condensation. Model-ing attempts have been generally unsuccessful in generatingthe atmospheric extensions of red supergiants necessary toput enough material at the dust sublimation front to re-cover the observed mass-loss rates of red supergiants (e.g.,Arroyo-Torres et al. 2015).An alternative suggestion has been that pulsational mo-tions might be accompanied or replaced by significant atmo-spheric turbulence (e.g., Gustafsson & Plez 1992; Josselin& Plez 2007), and that this turbulence might be seeded bythe vigorous convection expected in the atmospheres of redsupergiants (see Freytag et al. 2012, for a theoretical treat-ment of convection in red supergiant envelopes). At themoment, these convection simulations produce only mod-estly extended atmospheres compared to what is inferredfrom observations (e.g., Arroyo-Torres et al. 2015). Never- Article number, page 1 of 14 a r X i v : . [ a s t r o - ph . S R ] J a n &A proofs: manuscript no. rsg_modParker_final theless, observations of red supergiants do suggest that theouter 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 thederivation of an analytic, theoretical mass-loss rate that fo-cuses on these observed turbulent velocities present in theatmospheres of red supergiants. Doing so, we find that theinferred turbulent motions are, alone, sufficient to explainthe mass-loss rates of red supergiants even in absence of anydust opacity. This model then can be seen as a turbulentpressure driven extension of the classical thermal pressuredriven Parker wind (Parker 1958). In contrast to the so-lar wind and other applications of this theory (for instancethe “warm wind” model, Hearn 1975; Rogerson & Lamers1975), however, the total sound and turbulent speed can bemuch lower due to the lower surface gravity of red super-giants. In Sect. 2 we lay out the basic formalism for thismodel in an isothermal set-up, under which the mass-lossrate is fully analytic. By imposing a temperature structure,we then extend the model in Sect. 3 by iteratively solvingthe relevent equations to converge to non-isothermal mass-loss rates. We further fit over the difference between theisothermal and non-isothermal mass-loss rates to provide atheoretically motivated mass loss rate as an analytic func-tion of stellar parameters and turbulent velocity. In Sect. 4we review some of the currently used empirical mass-lossrate prescriptions for red supergiants and compare our the-oretical rate to these, as well as to the sample of observedred supergiants from Josselin & Plez (2007) and the ob-servations of Antares by Ohnaka et al. (2017). Finally, wediscuss some future directions for this model in Sect. 5.
2. Analytic mass-loss rate from levitation of anisothermal atmosphere
To begin, we first consider the relevant equations that mustbe satisfied at all points in an isothermal flow, namelyconservation of mass and momentum. For a purely radial,spherically symmetric outflow in a steady state these are v ∂ρ∂r + ρ ∂v∂r + 2 ρvr = 0 (1) v ∂v∂r = − ρ ∂P∂r − GM ∗ r , (2)expressed in terms of radial velocity v , density ρ , and totalpressure P as functions of spherical radius r , as well asstellar mass M ∗ and gravitational constant G .As is done in standard static model atmosphere com-putations (see, e.g., Gustafsson et al. 2008; their Sect. 3)pressure is expressed as the sum of thermal pressure P therm = ρc , turbulent pressure P turb = ρv , and ra-diation pressure P rad where the characteristic velocities ofthe first two are respectively the sound speed c s and theturbulent speed v turb . As concerns the radial componentof radiation pressure gradient, we associate this with theacceleration g rad such that − ρ ∂P rad ∂r = g rad = κF r c = κL ∗ πr c = Γ GM ∗ r , (3) with c the speed of light, κ the flux weighted mean opac-ity [cm / g] , and F r the radial component of the flux, re-lated to the stellar luminosity L ∗ as F r = L ∗ / (4 πr ) .The final equality further introduces the Eddington factor Γ ≡ κL ∗ / (4 πGM ∗ c ) . Explicitly replacing pressure with thiscombination and substituting Eq. 1 into Eq. 2 to eliminatedensity gradient terms, we find v (cid:18) − c + v v (cid:19) ∂v∂r = 2 (cid:0) c + v (cid:1) r − GM ∗ (1 − Γ) r , (4)where we have assumed that v turb is constant. We discussthis approximation further below.We note that Eq. 4 is the same equation as is solved foran isothermal, pressure driven Parker wind (Parker 1958),only now with a modified “effective sound speed”, c s , eff ≡ (cid:112) c + v . At the location where the flow velocity reachesthis effective sound speed the left side of Eq. 4 goes to zero,and therefore the right side must also go to zero. Solvingfor this criteria yields that the location of this critical pointin the wind outflow is given by a modified Parker radius, R p , mod = G M ∗ (1 − Γ)2 ( c + v ) . (5)This then suggests that the problem of wind accelerationhere can be envisioned to consist of two parts. First, a low-speed near “levitation” of material to this modified Parkerradius at which v = c s , eff , and, second, a subsequent accel-eration of the material to infinity.Below the modified Parker radius the contribution ofthe advection term ( v ∂v/∂r ) in Eq. 2 is almost negligibleso that the equation reduces to the standard equation forhydrostatic equilibrium, ρ ∂ρ∂r = − G M ∗ (1 − Γ)( c + v ) r . (6)Typically, the Eddington factor and, as alluded toabove, the turbulent speed may be expected to vary in ra-dius. However, under the assumption that this variation isonly mild, we may replace the exact expression Γ( r ) with Γ and, as already done in Eq. 4, v turb ( r ) with v turb . Theserepresentative constant values of Γ and v turb then in gen-eral would designate some appropriate spatial means (seealso the discussion in Gustafsson et al. 2008). Then Eq. 6can be integrated analytically from the stellar radius R ∗ toan arbitrary radius r ≤ R p , mod using the boundary valuedefinition ρ ∗ ≡ ρ ( R ∗ ) to give ρ ( r ) = ρ ∗ exp (cid:20) − R ∗ H (cid:18) − R ∗ r (cid:19)(cid:21) , (7)where we have made the substitution HR ∗ ≡ R ∗ ( c + v ) G M ∗ (1 − Γ) = 2 c + v v , eff = 12 R ∗ R p , mod . (8)Here v esc , eff ≡ (cid:112) GM ∗ (1 − Γ) /R ∗ denotes the escape speedfrom R ∗ with √ − Γ accounting for the reduction in effec-tive gravity due to radiative acceleration. Article number, page 2 of 14. D. Kee et al.: Analytic RSG mass-loss rates from turbulent pressure
In order to determine the value of ρ ∗ , we make use ofthe definition of the stellar radius in optical depth space, τ ( R ∗ ) ≡ / , from a spherical extension of the classicalplanar gray model atmosphere (Lucy 1971; and see furtherSect. 3). Within this model, one defines a "spherically mod-ified" optical depth scale as τ ( r ) ≡ (cid:90) ∞ r κ ρ (cid:18) R ∗ r (cid:19) dr , (9)and, as such, κ (cid:90) ∞ R ∗ ρ ( r ) (cid:18) R ∗ r (cid:19) dr = 23 , (10)where we have again used the above approximation that Γ , and therefore κ , is independent of radius. Inserting thehydrostatic density profile from Eq. 7 into Eq. 9 yields ananalytic value of ρ ∗ , ρ ∗ = 23 κ H (cid:20) − exp (cid:18) − R ∗ H (cid:19)(cid:21) − . (11)Before combining the above terms to derive a turbulentpressure driven mass-loss rate, it is important to examine afew of the assumptions that went into the derivation. First,we consider the approximation that the advection term isnegligible for r < R p , mod . If we relax this approximation, wefind that integrating Eq. 2 from R ∗ to an arbitrary radius r gives log (cid:18) ρ ( r ) ρ ∗ (cid:19) = − R ∗ H (cid:18) − R ∗ r (cid:19) −
12 ( c + v ) (cid:90) rR ∗ ∂v ∂r dr . (12)While in general solving this equation at an arbitrary ra-dius requires numerical integration, the integration from R ∗ to R p , mod remains analytic as v ( R ∗ ) / ( c + v ) ≈ and v ( R p , mod ) / ( c + v ) ≡ . Comparing ρ ( R p , mod ) com-puted with and without the advection term shows that in-cluding this term reduces the density at the modified Parkerradius by a constant factor exp( − / . As this reduction isanalytic and constant, we therefore include it in all calcu-lations for the remainder of this section.Related to this change in density, it is also important toexamine whether the optical depth, and as such the basedensity ρ ∗ , is significantly impacted by the inclusion of theadvection term or by a radial profile in the opacity. To testthis we plot τ ( r/R ∗ ) computed from the hydrostatic den-sity structure with a constant opacity in Figure 1, assumingthat H/R ∗ = 0 . , a characteristic value for a red super-giant atmosphere. We note here, from examining Eqns. 7, 9,and 11, that τ ( r/R ∗ ) only depends on this chosen H/R ∗ asthe definition of τ (1) = 2 / results in κ in Equations 9 and11 cancelling. As expected, almost all of the optical depthis accumulated within only a few scale heights of the stellarsurface, here within about half a stellar radius of R ∗ . This iswell away from the modified Parker radius, here ∼ R ∗ , andtherefore the optical depth used to define the stellar radiusand ρ ∗ is all accumulated in a region where the advectionterm is indeed negligible. Similarly, this means that κ and r ( R )10 Fig. 1: Optical depth calculated inward to r. Note that mostof the optical depth is accumulated within a stellar radiusof the star. ( v + c ) only need to be constant over this same smallspatial extent in order for ρ ∗ to be unaffected by any vari-ations they may have. Sect. 3 will revisit the implicationsof allowing c s to take on a realistic radial profile.Finally, to turn the terms derived thus far into a mass-loss rate, we note that Eq. 1 is equivalent to the constraint ρvr = const.; the total mass-loss rate is this conservedquantity integrated over all solid angle, ˙ M ≡ πρvr . Fur-ther, as this quantity is independent of radius, we can sim-ply take its value at the modified Parker radius ˙ M = 4 π ρ ( R p , mod ) (cid:113) c + v R , mod , (13)as the total mass-loss rate of the star. Combining Eqs. 7, 8,and 11 and accounting for the additional factor exp( − / from advection gives ρ ( R p , mod ) to be ρ ( R p , mod ) = 43 R p , mod κR ∗ exp (cid:104) − R p , mod R ∗ + (cid:105) − exp (cid:104) − R p , mod R ∗ (cid:105) . (14)All quantities in Eqs. 13 and 14 are analytically knownwithin this formalism, such that these expressions togetheroffer a fully analytic mass-loss rate. In understanding the regimes of stellar parameters wheresuch a mass-loss rate model is and is not viable, it is im-portant to note that
H/R ∗ = R ∗ / (2 R p , mod ) as shown inEq. 8. This means that the same scale which controls theexponential stratification of density also controls where thecritical point lies. Therefore, H/R ∗ must be relatively largecompared to, for instance, the Sun in order to prevent theexponential density stratification of the atmosphere fromdriving ˙ M to zero. At the same time, R p , mod and by exten-sion R ∗ must be large (in absolute units) simply due to the R , mod dependence of ˙ M . Article number, page 3 of 14 &A proofs: manuscript no. rsg_modParker_final M ( M )10 M ( M / y r )
500 750 1000 1250 R ( R )10 M ( M / y r ) T eff (K)10 M ( M / y r ) Fig. 2: Dependence of the analytic mass-loss rate on stellar mass (left), stellar radius (center), and temperature (right).The mass-loss rate derived in Sect. 2.1 is a function ofonly a few key parameters of the physical set-up, namelystellar mass M ∗ , stellar radius R ∗ , sound speed c s , turbulentvelocity v turb , and opacity κ . Specifically, ˙ M ∝ κ (cid:18) R ∗ H (cid:19) exp (cid:18) − R ∗ H (cid:19) R ∗ , (15)where H/R ∗ itself scales as HR ∗ ∝ R ∗ (cid:0) c + v (cid:1) M ∗ (1 − Γ) . (16)We can trade the sound speed dependence for a stellarparameter by assuming for the following discussion thatthe wind has c s = (cid:112) k B T eff /m H with Boltzmann constant k B , mean molecular mass equal to the mass of a hydro-gen atom m H , and temperature equal to the stellar effec-tive temperature T eff = ( L/ (4 π R ∗ σ )) / , where σ is theStefan-Boltzmann constant. To investigate the behavior ofthis analytic expression, we begin by calculating the mass-loss from a star with parameters selected to be consistentwith a “typical” RSG, namely M ∗ = 10 M (cid:12) , R ∗ = 500 R (cid:12) , T eff = 3500 K ( c s = 5 . km s − ), v turb = 15 km s − , and κ = 0 . cm g − . Note that here opacity is estimated fromthe OPAL opacity tables (Iglesias & Rogers 1996) as a lowerlimit to the total opacity as will be discussed further be-low. For these parameters, ˙ M = 8 . × − M (cid:12) yr − , inthe range of plausible mass-loss rates for a red supergiant(see, e.g., Josselin & Plez 2007).For the remaining discussion in this section, we vary in-dividual model parameters while holding all others to thevalues provided immediately above to show how ˙ M reactsto such variations. The scaling of this mass-loss rate withstellar parameters is relatively straightforward as shown byFigure 2. Specifically, increasing the stellar radius or effec-tive temperature over a reasonable range for red supergiantsincreases the mass-loss rate because increasing either ofthese increases H/R ∗ . This goes into the exponential termin Eq. 15 which increases faster than the power-law termsfall off. Conversely, increasing M ∗ has the opposite effect on H/R ∗ and therefore increasing stellar mass decreases mass-loss rate. Note that T eff only appears through the soundspeed, and thus only in the combination c + v . There-fore, it is unsurprising that varying this parameter has amuch weaker effect than either stellar mass or radius, afact we return to later in the discussion.Next, we examine the scaling of mass-loss rate with v turb , shown in Figure 3. In order to emphasize the im-portance of v turb as the driving mechanism for this model, v turb (km/s)10 M ( M / y r ) Fig. 3: Dependence of the analytic mass-loss rate on tur-bulent velocity.we vary v turb from a minimum microturbulence for red su-pergiants ( ∼ − ) up through the characteristic tur-bulent velocities inferred by Josselin & Plez (2007). As wasthe case with the variation of stellar parameters, for muchof the range of v turb the dominant effect is on the scaleheight. As H/R ∝ v for v turb (cid:29) c s , this dependenceis quite steep. However, note here that we allow a muchlarger range of variation in v turb than in the other parame-ters, which means that we see an additional regime that didnot appear in the variation of stellar parameters, namelythe asymptotic mass-loss rate reached at high v turb . Here exp( − R ∗ /H ) begins to vary less rapidly with the increasein H/R ∗ and the decrease of R p , mod and ρ ∗ with increas-ing H/R ∗ takes over. The scaling is truncated at the pointwhere increasing v turb means R p , mod < R ∗ , as the model isno longer meaningful in this regime. If v turb were allowedto vary to arbitrarily low values, a lower asymptote wouldalso appear corresponding to the negligibly low mass lossthat could be driven by c s alone.Finally, we examine the scaling of ˙ M with opacity,shown in Figure 4. While the first calculation of ˙ M forcharacteristic red supergiant parameters above assumed theRosseland mean opacity from OPAL, the acceleration of thewind means that a purely static opacity model may not bethe most accurate one (see also Gustafsson & Plez 1992).In fact, the Doppler shifting of spectral lines from atomicand molecular species will increase the characteristic opac- Article number, page 4 of 14. D. Kee et al.: Analytic RSG mass-loss rates from turbulent pressure M ( M / y r ) M ( M / y r ) M ( M / y r ) Fig. 4: Dependence of the analytic mass-loss rate on Γ for M ∗ = 10 M (cid:12) and T eff = 3500 K. The left panel uses R ∗ = 500 R (cid:12) and v turb = 10 km s − , the middle panel uses R ∗ = 500 R (cid:12) and v turb = 15 km s − , and the right panel uses R ∗ = 1000 R (cid:12) and v turb = 15 km s − .ity scale of the wind, as discussed in Appendix A. In prac-tice, the computation presented in Appendix A for κ and Γ from these Doppler shifted spectral lines becomes computa-tionally expensive to include in a model such as the one wediscuss in the remainder of the paper. Specifically, millionsof spectral lines are involved (see especially the tabulationsof infrared spectral lines of water, e.g., Barber et al. 2006),with the opacity of each line becoming position-dependentdue both to the velocity gradient and the radial temper-ature variation. In studies of hot, OB star winds, this iscircumvented by introducing a distribution function thatallows one to analytically approximate the cumulative ef-fects of all these spectral lines. As a first step, we have takenthe brute force approach of computing the full sum but onlyfor fixed hydrodynamic structures, thus allowing us to ap-proximate the net effect of this additional radiative driving.By computing the force from infrared spectral line lists ofCO (Goorvitch 1994), TiO (McKemmish et al. 2019), andH O (Barber et al. 2006) as a sample, we find that the netincrease in flux-weighted mean opacity may be an order ofmagnitude or more over the basal OPAL opacity. However,this preliminary study deserves further attention as it isbased on analyses of fixed hydrodynamic structures, whileline acceleration is a sensitive function of the velocity field.Further, we have only included three molecules here, andalso used a simplified radiative transfer treatment based onthe so-called Sobolev approximation (see Appendix A), im-plying the real effect could be a bit different than the 1 dexmentioned above.Here we take a more general view of the effect ofopacity variation, simply allowing Γ to increase from itsbasal value as implied by the Rosseland mean opacity, Γ Ros . As was the case for v turb , we truncate the plotwhen R p , mod = R ∗ as this is again the limit at whichthe model breaks down. To emphasize the complex be-havior induced by varying Γ and thus κ , we show thescaling of ˙ M for three models: v turb = 10 km s − and R ∗ = 500 R (cid:12) (left panel of Figure 4), v turb = 15 km s − and R ∗ = 500 R (cid:12) (center panel of Figure 4), v turb = 15 km s − and R ∗ = 1000 R (cid:12) (right panel of Figure 4). Note here that Γ Ros = κ R L ∗ / (4 π G M ∗ c ) = κ R R ∗ σ T / ( G M ∗ c ) . All pa-rameters are the same for the three cases except R ∗ , suchthat this minimal point for Γ is lower for the two modelswith R ∗ = 500 R (cid:12) ( Γ Ros = 2 . × − ) than the one with R ∗ = 1000 R (cid:12) ( Γ Ros = 1 . × − ).In general, as increasing Γ means increasing the radia-tion force available to levitate the atmosphere, one might expect that ˙ M would monotonically increase with Γ . How-ever, it is important to recall that we choose our uniquesolution for the mass-loss rate by requiring τ ( R ∗ ) = 2 / through Equation 10, which results in an inverse depen-dence of the stellar photosphere density on opacity in Equa-tion 11. When Γ is small and the wind is predominantlydriven by turbulent pressure, this inverse dependence winsout and increasing Γ and therefore κ implies a reduction ofthe density of the wind and therefore also of the mass-lossrate to keep R ∗ fixed.As Γ approaches unity, radiation then contributes moremeaningfully to the total force budget of the wind and theexpected increase manifests itself. In the case that the tur-bulent pressure gradient is not already driving a strongwind mass loss (e.g. v turb and R ∗ small as in the left panelof Figure 4) then radiative acceleration is able to drive amuch stronger mass loss near Γ unity than pressure could atlow Γ . In the case where turbulent pressure already drivesa strong mass loss (i.e., v turb and R ∗ large as in the rightpanel of Figure 4) the net effect is only to cancel out thereduction in ˙ M from the increased opacity. However, eventhis increase (or flattening) is not able to hold out all theway to Γ unity. Eventually, as was the case with increas-ing v turb , R p , mod is driven to R ∗ and ˙ M drops off. Again,this point merely denotes where such a turbulent pressuredriven wind model breaks down.As mentioned in the introduction, it is intriguing thatthis model is able to recover reasonable mass-loss rates forred supergiants without appealing to dust opacity. In fact,as Γ Ros is small in the standard case we consider, this modeldoes not appeal to significant radiative driving from anysource. One of the theories present in the literature is thatred supergiant winds behave analagously to the winds ofAGB stars with some atmospheric extension allowing forthe condensation of dust, the opacity of which is the mainmass-loss rate driver (see e.g. Höfner & Olofsson 2018, fora review). Instead, we here find that it is plausible thatthe extreme extension of their atmospheres by turbulentpressure alone could be the dominant part determining thered supergiant mass-loss process.However, it is important to note here that opacity maystill play an important role in altering the structure of thewinds of red supergiants. While we have here assumed aconstant opacity as a function of distance from the star inorder to examine the scaling of ˙ M with Γ , this is not ex-pected to be the case. As the wind cools away from the Article number, page 5 of 14 &A proofs: manuscript no. rsg_modParker_final star, additional molecules and dust will form, thereby en-hancing the wind opacity. Even the continuum opacity islikely to shift as the Hydrogen anion is no longer the dom-inant continuum opacity source. As discussed in the priorsubsection, the effect of this enhanced opacity on at leastthe optical depth scale, and as such the density scale usedin this model, may not be substantial as long as the opac-ity enhancement occurs beyond the first stellar radius orso. This is because, as shown by Fig. 1, the majority ofoptical depth is accumulated in a region quite near to thestar making κ in Eqn. 11 effectively an average over thisnear star region. Enhancing the opacity further out would,however, contribute to extending the scale height of theatmosphere, and even potentially allow the wind to switchfrom being turbulence driven to radiation driven at some lo-cation. Both these effects could increase ˙ M somewhat overwhat is presented here. In light of these open items relatedto wind opacity, as we discuss further in Sect. 5, futurework should re-examine what role depth dependent dust,molecular, and continuum opacity may play, especially insetting the terminal velocity of the stellar wind as this is ill-defined under a Parker-like pressure driven wind (see alsothe discussion in Sect. 3). For the remainder of this paper,however, given the highly promising results throughout, weproceed with turbulent pressure supplemented only by themarginal additional contribution of a constant Γ Ros .
3. Numerically determined mass-loss rate fromsteady-state, nonisothermal, dynamical winds
While the simplified, isothermal treatment has the ben-efit of being purely analytic, we can also treat the non-isothermal wind structure in a numerical approach. To doso, we use the spherically symmetric generalization to theplane parallel gray atmosphere introduced by Lucy (1971) T = T (cid:32) W + 34 (cid:90) rR ∗ κ ρ (cid:18) R ∗ r (cid:19) dr (cid:33) (17) = T (cid:18) W + 34 τ (cid:19) , (18)where W ≡ (1 − (cid:112) − ( R ∗ /r ) ) / is the geometric dilutionfactor. The second equality again uses the spherically mod-ified optical depth as given by Eq. 9. For low optical depthmaterial, and for distances away from the stellar surface thistemperature structure behaves like T ∝ / √ r as is gener-ally assumed for the dust free regions near red supergiants(e.g. Decin et al. 2006). Accounting for this temperaturestructure in Equations 1 and 2 while still maintaining theassumption that v turb is constant leads to the system ofnon-isothermal equations v ∂ρ∂r + ρ ∂v∂r + 2 ρ vr = 0 (19) v ∂v∂r = − c + v ρ ∂ρ∂r − G M ∗ (1 − Γ) r − k B m H ∂T∂r (20) ∂T ∂r = T (cid:32) ∂W∂r + 34 κρ (cid:18) R ∗ r (cid:19) (cid:33) . (21) Note that, by combining this system of differential equa-tions into a single equation to eliminate ∂ρ/∂r and ∂T /∂r in Equation 20, we can see that these equations still containthe same critical point where v = c + v , although nowwith c s a function of r .In order to solve these differential equations, we beginby assuming initial velocity and density profiles v ( r ) = (cid:112) c + v − R p , mod /r (cid:18) − R ∗ r (cid:19) (22) ρ ( r ) = ˙ M π v ( r ) r . (23)As boundary conditions, we fix R ∗ and impose that T ( R ∗ ) = T eff . These are then sufficient conditions to be-gin the following iteration procedure.1. A radial grid is built up inward (toward R ∗ ) and out-ward (away from R ∗ ) from the critical point, R p , mod , asdetermined in the prior iteration (or the initial condi-tions for the first iteration). Radial points are spaced by H ( r ) / to properly resolve the variations in all hydro-dynamic quantities.2. The prior structures in velocity, density, and tempera-ture are interpolated onto the new radial grid.3. Gradients of all necessary quantities for Eqns. 19, 20,and 21 are computed at each grid point. At the criticalpoint, r = R p , mod , where the velocity gradient is ill-defined, l’Hopital’s rule is used.4. Using these gradients, Eqns. 19, 20, and 21 are numeri-cally solved at each point on the radial grid inward andoutward from the critical point, R p , mod .5. The resulting velocity, density, and temperature as func-tions of radius are used to begin the process again fromStep 1.This loop is repeated until velocity, density, and tempera-ture at each point agree from one iteration to the next tobetter than 1%. At this point, the total optical depth ofthe wind is computed at the current mass loss rate. As thestellar radius is defined to be at τ ( R ∗ ) = 2 / , the stellarmass-loss rate is re-computed according to ˙ M new = ˙ M old / τ old ( R ∗ ) . (24)As optical depth is proportional to density, which is in turnproportional to ˙ M , a constant change in mass-loss rateand thus density at every point has the effect of forcing τ ( R ∗ ) = 2 / , assuming v ( r ) is unchanged. As changing ρ ( r ) will certainly change v ( r ) , however, the final velocity pro-file is combined with the updated mass loss as new initialconditions to restart the loop described above. This outeriteration to update mass-loss rate is repeated until the up-date in ˙ M is also less than 1%, at which the final mass lossrate is returned as ˙ M num . Note that for the remainder ofthe paper we include the subscript “num” to differentiatethis numerically determined, non-isothermal mass-loss ratefrom the isothermal, analytic one derived in Sect. 2, whichwe henceforth call ˙ M an .Here it is again interesting to examine our “typical”model with M ∗ = 10 M (cid:12) , R ∗ = 500 R (cid:12) , T eff = 3500 K, Article number, page 6 of 14. D. Kee et al.: Analytic RSG mass-loss rates from turbulent pressure R )10 ( g / c m ) R )0102030 v ( k m / s ) R )10 T ( K ) Fig. 5: Radial dependence of density (left), velocity (center), and temperature (right) for the non-isothermal case with M ∗ = 10 M (cid:12) , R ∗ = 500 R (cid:12) , T eff = 3500 K, v turb = 15 km s − , and κ = 0 . cm g − . The red dashed line denotes themodified Parker radius.
14 16 18 20 22 24 v turb (km s )0.20.30.40.50.60.70.80.9 M nu m / M a n M ( M ) Fig. 6: Ratio of ˙ M calculated numerically for a non-isothermal situation to the analytic ˙ M calculated for anisothermal wind plotted as a function of v turb and coloredby M ∗ , varying R ∗ , M ∗ , and v turn . The effective tempera-ture and Rosseland opacity are held fixed at T eff = 3500 Kand κ R = 0 . cm g − . v turb = 15 km s − , and κ = 0 . cm g − . Using the methoddescribed above, we find that it has a non-isothermal mass-loss rate of ˙ M num = 4 . × − M (cid:12) yr − , approximatelyhalf the analytic rate. We can also take this opportunity toplot the radial profiles of density, velocity, and temperatureas shown in Figure 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 predominantlyaccumulated below R p , mod such that the mass-loss rate ineach simulation is primarily determined by integration overonly these near star regions. However, like the isothermal,thermal pressure-driven Parker wind, the total driving en-ergy available to the wind diverges as the outer radius goesto infinity because v turb is constant, such that v ∞ is notuniquely determined by such a model. While the terminalvelocity of the wind will depend on the details of the radial v turb profile, radial temperature profile, and any additionalforces acting on the outer wind (for instance dust opacity),it is promising that these simulations already generate ve-locities at 30 R ∗ comparable to typically inferred values of v ∞ for red supergiants. While we here present the mass-loss rate from a simpli-fied gas temperature structure arising from radiative equi-librium, in reality red supergiant winds may have regionsthat are heated by shocks and dissipation of turbulence(see, e.g. Schrijver & Zwaan 2000). In general, such shockheating and dissipation is not inferred to generate windtemperatures throughout the full wind acceleration regionsignificantly in excess of the stellar effective temperature,nor should the wind be expected in general to cool be-low radiative equilibrium. Therefore, the simplified cases wepresent here of an isothermal wind at the stellar effectivetemperature and a wind with temperature set according tothe spherically modified gray atmosphere can be roughlytaken as thermal structures leading to, respectively, upperand lower limits to the mass-loss rate. For individual ob-jects where more is known about the temperature profile,however, the method presented in this section can be di-rectly applied by replacing the temperature information inEquations 17 and 21. As a sample case, let us take the chro-mospheric temperature profile inferred by O’Gorman et al.(2020) for Antares. In these observations, gas temperatureis inferred to decrease from the stellar effective tempera-ture to a minimum around . R ∗ before coming back upto a maximum approximately equal to the stellar effectivetemperature around . R ∗ , and then finally falling off out-ward from there (see O’Gorman et al. 2020, their Figure3). Inserting this temperature profile into our model withall other parameters still as defined in our “typical” modelabove returns a mass-loss rate . × − M (cid:12) yr − . As antic-ipated, this falls between the limiting mass-loss rates froman isothermal model at the stellar effective temperature( . × − M (cid:12) yr − ) and a model using the Lucy (1971)spherically modified gray atmosphere ( . × − M (cid:12) yr − ). To investigate the role of the radial decline in temperatureon the mass-loss rate, we run a grid of steady-state modelsvarying stellar mass, stellar radius, and turbulent velocityand then plot the ratio ˙ M num / ˙ M an . For this test we fix T eff = 3500 K, motivated by the only modest dependenceof ˙ M on this property (see Figure 2), and continue to use κ R = 0 . cm g − . As shown by Figure 6, this correctionfactor is itself a function of the model parameters that wecan fit. We choose two methods of doing so, one where wefit the product of a constant normalization and power-laws Article number, page 7 of 14 &A proofs: manuscript no. rsg_modParker_final
Table 1: Comparison of the fit methodsMethod 1 Method 2 χ (cid:104) Error (cid:105) σ Error v turb /v esc ( M ∗ , R ∗ ) (hereafter Method 2).As figures of merit to compare the two fits, we use χ ,the mean and standard deviation of the difference betweenthe fit formula and the actual ratio ˙ M num / ˙ M an , and themaximum error between the fit and ˙ M num / ˙ M an . The com-parison of the best models from each method is summarizedin Table 1. Given that method 1 uses four free parameters(a mean and three power-law indices) versus the two freeparameters used by method 2 (a mean and one power-lawindex), and that the figures of merit we have selected tocompare between the two methods are comparable, we pro-pose that the two parameter fit from method 2 should bepreferable moving forward. Therefore, the non-isothermalcorrection factor is (cid:32) ˙ M num ˙ M an (cid:33) = (cid:18) v turb / (17 km s − ) v esc ( M ∗ , R ∗ ) / (60 km s − ) (cid:19) . , (25)where we have omitted the leading constant as the best fitfor this constant was 10 . = 1 . .Note from Table 1 that this fit is quite good, with thecorrection factor recovered to about a factor of two in theworst case for the two parameter fit. Therefore we presentthe combination of Equation 25 with Equations 5, 13, and14 as a predictive, theoretical mass-loss rate scaling for redsupergiants.
4. Comparison of the theoretical mass-loss ratescaling with empirical rates and observedmass-loss
Now that we have this theoretical, predictive mass-loss rateas a function of stellar parameters, a natural next step is tocompare this with empirical fits to red supergiant mass lossand individual observations. Given that no systematic the-oretical predictions for RSG mass loss exist at the present,such empirical fits constitute the state of the art regardingRSG mass loss in applications such as stellar evolution. Wereview the empirical rates implemented in the stellar evo-lution code MESA, as well as the recently published ratesfrom Beasor et al. (2020) in Sect. 4.1 in order to build up arepresentative sample. In Sect. 4.2 we then compare theseempirical rates to our theoretical ones. Here we highlightthat the empirical rates in general do not depend on allthree fundamental stellar parameters (i.e. M ∗ , R ∗ , and L ∗ )while the theoretical mass-loss rates presented here does.Therefore Sect. 4.2 also includes a discussion of the scal-ings of stellar luminosity, mass, and radius with one an-other that we use to perform the comparison of our rates with the empirical ones. Finally, in Sect. 4.3 we compareour predicted rates with observationally inferred mass-lossrates for a selection of red supergiants. The current version of the stellar evolution code
MESA (Pax-ton et al. 2011, 2013, 2015, 2018, 2019) provides built-in ac-cess to three empirical mass-loss rates for red supergiants.These come from de Jager et al. (1988) ˙ M = 10 − . (cid:18) L ∗ L (cid:12) (cid:19) . ( T eff ) − . , (26)Nieuwenhuijzen & de Jager (1990) ˙ M = 10 − . (cid:18) L ∗ L (cid:12) (cid:19) . (cid:18) M ∗ M (cid:12) (cid:19) . (cid:18) R ∗ R (cid:12) (cid:19) . (27)and van Loon et al. (2005) ˙ M = 10 − . (cid:18) L ∗ L (cid:12) (cid:19) . (cid:18) T eff (cid:19) − . . (28)To these we add the recently published empirical rate fromBeasor et al. (2020) ˙ M = 10 − . − . M ∗ /M (cid:12) (cid:18) L ∗ L (cid:12) (cid:19) . , (29)as the authors identify that this rate is notably lower thanthose previously published, a fact that we use to provide anindication to the range of possible observationally inferredrates. While these four are far from the only empirical mass-loss prescriptions that are used, they provide a representa-tive sample. For a review of a variety of other empiricalrates, see, for example, Mauron & Josselin (2011). As mentioned in the prior subsection, it is common for em-pirical rates to be provided as functions that do not in-clude a dependency on all three fundamental stellar pa-rameters (i.e. M ∗ , R ∗ , and L ∗ , or equivalently M ∗ , L ∗ , and T eff ). For the theoretical rates, however, all three parame-ters are needed. To obtain these, we consider stellar evo-lution tracks for nonrotating stars from both MESA (Pax-ton et al. 2011) and
GENEC (Ekström et al. 2012). Notethat
GENEC uses a somewhat different prescription for themass-loss rate than what is discussed above for
MESA . For M ∗ ≤ M (cid:12) GENEC adopts the Reimers (1975, 1977) rates,whereas for M ∗ ≥ M (cid:12) GENEC uses de Jager et al. (1988)for log ( T eff ) > . and a linear fit to data from Sylvesteret al. (1998) and van Loon et al. (1999) otherwise. Despitethese differences in mass-loss prescription, comparison ofPaxton et al. (2011) and Ekström et al. (2012) shows thatin both cases the red supergiant branch roughly follows L ∗ = f M ∗ as shown by Figure 7 with f = 12 . for MESAand f = 18 . for Geneva. Therefore, as a compromise wechose f = 15 . for our mass-luminosity relation to comparewith the empirical rates. As a simplifying assumption, we Article number, page 8 of 14. D. Kee et al.: Analytic RSG mass-loss rates from turbulent pressure
107 8 9 20 M ( M )10 L ( L ) GenevaMESA
Fig. 7: Luminosity as a function of mass from red super-giants in Paxton et al. (2011) and Ekström et al. (2012).The overplotted curves assume L ∗ = f M ∗ with f = 12 . for the MESA data and f = 18 . for the Geneva data.continue to hold T eff fixed at 3500 K as was done in Sect. 3,so that L ∗ ∝ R ∗ .Using these mass-radius-luminosity relationships, wecan now compare the empirical and theoretical rates asshown by Figure 8. As the empirical rates differ from oneanother by several orders of magnitude at all luminosities,we select four values of v turb from 15 to 21 km s − for thiscomparison. Despite this scatter between individual empiri-cal rates and in turn between the theoretical rates using dif-ferent turbulent velocities, the theoretical rates are broadlyconsistent 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 em-pirical 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 su-pergiants. Alternatively, as is suggested by the mismatchin both value and slope of the empirical rates with one an-other, this shallower slope may simply reflect some missingphysics when inferring and calibrating these empirical ratesor potential biases in the observed sample.We note that the mass-luminosity-radius relations thatwe use here may not be appropriate to Galactic stars start-ing out their lives rapidly rotating, as these tend to evolveto 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 temper-atures of these stars on the red supergiant branch this isconsistent with these stars evolving to have larger radii,lower surface gravities, and higher mass-loss rates if allother physics remains the same. This further suggests thattwo stars with different masses can reach the same point onthe HR diagram with different mass-loss rates which wouldconfound any mass-loss rate predictions omitting a depen-dence on M ∗ . In fact the general treatment that mass-lossrates are enhanced by rotation would go in the direction ofonly 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 mass-loss rate with increased L ( L )10 M ( M / y r ) de JagerNieuwenhuijzenvan LoonBeasorTheoretical, v turb = 15 km/sTheoretical, v turb = 17 km/sTheoretical, v turb = 19 km/sTheoretical, v turb = 21 km/s Fig. 8: An overplot of three of the theoretical mass-lossrate scalings used in MESA with the predicted mass-lossrate from the theoretical scaling. Four values of v turb areused as noted in the plot label.rotation rate are based on fits to the theoretical predictionsof Friend & Abbott (1986). This work addresses only thewind mass-loss rate of hot stars with radiation driven out-flows, and has been questioned by subsequent work (e.g.Owocki et al. 1996; Petrenz & Puls 2000) such that thequestion of the role of rotation in the mass-loss of red su-pergiants remains open. While the preceding comparison of the empirical rates withtheoretical ones helps to generally motivate the quality ofthe theoretical mass-loss rate predictions, we can go onestep further and compare theory and observation for a sam-ple of observed stars. For this we select the observationscompiled by Josselin & Plez (2007) with the addition ofAntares using data from Ohnaka et al. (2017). Note thatthis limited sample is based on these authors compiling ve-locity dispersion measurements that we can compare withthe v turb required by our model. Table 2 recapitulates thestellar masses, effective temperatures, and radii from Jos-selin & Plez (2007) and Ohnaka et al. (2017). As the modelwe present predicts gas mass-loss rates while Josselin & Plez(2007) report dust mass-loss rates, we apply a constant gasto dust mass ratio ˙ M gas = 250 ˙ M dust . This ratio is itself notparticularly well constrained by the existing literature asshown by a sample comparison of the gas mass-loss rates ofDe Beck et al. (2010) with the dust mass-loss rates of Har-wit et al. (2001) and Verhoelst et al. (2009) for VY CMa, µ Cep, and α Ori, which returns ratios between ∼
70 and ∼ v turb , Obs , and the tur-bulent velocity that would be required for the theoreticalmodel to reproduce ˙ M gas as v turb , Theory . For the data fromJosselin & Plez (2007), we round their reported values tothe nearest km s − as the errors on the original data maybe as large as ± km s − . We report our required v turb withthe same accuracy. Article number, page 9 of 14 &A proofs: manuscript no. rsg_modParker_final
Table 2: Table of observed parameters from Josselin & Plez (2007) with the addition of Antaresas observed by Ohnaka et al. (2013, 2017), including v turb as required by theory to recover themass-loss rate keeping all other parameters fixed.Number Name Mass T eff Radius ˙ M gas a v turb , Obs v turb , Theory M (cid:12) K R (cid:12) − M (cid:12) yr − km s − α Ori 15 3780 589 5.0 19 172 V466 Cas 12 3780 331 0.5 12 193 AD Per 12 3720 457 2.0 21 174 FZ Per 12 3920 324 1.75 16 205 BD+243902 15 4240 427 7.25 23 216 BI Cyg 20 3720 851 10.25 23 167 BC Cyg 20 3570 1230 8.0 22 138 RW Cyg 20 3920 676 8.25 20 199 SW Cep 9 3570 234 11.5 24 2310 µ Cep 25 3750 1259 3.75 23 1411 ST Cep 9 4200 174 6.25 23 2612 TZ Cas 15 3670 646 9.5 17 1713 Antares 12.7 b b b c d a A constant gas to dust mass ratio 250 has been assumed to convert ˙ M dust as compiled byJosselin & Plez (2007) to ˙ M gas . b Ohnaka et al. (2013) c Braun et al. (2012) d Ohnaka et al. (2017)
12 14 16 18 20 22 24 26 v (km s )012345 c o un t v turb, Obs v turb, Theory v t u r b ( k m / s ) v turb, Obs v turb, Theory Fig. 9: Comparisons of v turb from observations ( v turb , Obs ) and from matching corresponding observationally inferredmass-loss rates to the theoretical predictions presented here (see text). The left panel presents a histogram of thedistribution while the right panel plots v turb , Obs and v turb , Theory for each object numbered according to Table 2.Comparing the mean value of the observed turbu-lent velocity v turb , Obs to the mean theoretical turbu-lent velocity v turb , Theory , the two are consistent witheach other over the population to within one stan-dard deviation, as v turb , Obs = 20 . ± . km s − and v turb , Theory = 18 . ± . km s − . Repeating this ex-ercise with the extremal gas to dust mass ratios abovealso retrieves v turb values consistent with the v turb , Obs and v turb , Theory values reported here to within one standarddeviation, such that this choice has not significantly im-pacted our results. These ranges are also generally consis-tent with the turbulent velocities required to recover theempirical rates in the prior subsection. Moreover, the over-all ranges are consistent with v turb , Obs ∈ [12, 24] km s − and v turb , Theory ∈ [13, 26] km s − . Finally, looking at thedistributions of v turb for both theory and observation, plot-ted in the left panel of Figure 9, as well as comparing v turb , Obs to v turb , Theory for individual objects as plottedin the right panel of Figure 9, one can see that the major-ity of turbulent velocities required by theory (10 out of 13)are lower than what is observationally inferred. This arguesthat turbulent pressure is an ample driving source for thewind of red supergiants, strongly reinforcing the quality ofthis model in explaining the general behavior of these ob-jects’ mass loss.As turbulent velocity is not observationally inferred formany red supergiants, v turb = 18 . km s − as the mean v turb , Theory found above may thus be a reasonable value
Article number, page 10 of 14. D. Kee et al.: Analytic RSG mass-loss rates from turbulent pressure to use for application of our model in future works. Sim-ilarly, a reasonable range of variations on this would be14.8 km s − < v turb < − as suggested by plus orminus one standard deviation about this mean. While thiswill not provide a perfect match to the mass-loss rates ofindividual objects (see for example the right panel of Fig-ure 9) it provides a representative prediction for the averagemass-loss behavior of red supergiants such that it would beappropriate to use in stellar evolution codes. As such, andas discussed further in the following section, future workshould be dedicated to improving this estimation, estab-lishing whether v turb varies in time for individual stars, andtesting its dependence on key stellar parameters.
5. Summary and future directions
We have here derived a theoretical, analytic (dust-free)model for the mass loss-rates of red supergiants. Buildingupon the works by Gustafsson & Plez (1992) and Josselin& Plez (2007), we examine the role of turbulent pressure inlevitating material in the stellar atmosphere up to a mod-ified equivalent of the standard Parker wind critical point.This is achieved by employing a standard decomposition ofpressure into thermal and turbulent pressure componentsas, for instance, employed in standard 1D spherically sym-metric model atmospheres (e.g. Gustafsson et al. 2008).By applying the resulting expression for pressure in thesteady-state equations for conservation of mass and mo-mentum, we examine both hydrostatic and low velocitylevitation for an isothermal model (Sect. 2). We also nu-merically solve the full 1D steady-state equation-of-motionfor a non-isothermal model employing the Lucy (1971) tem-perature structure as a stand in for an energy equation. Bycombining each of these cases of steady-state levitation withthe constraint that the total accumulated optical depth ofcircumstellar material from infinity down to the stellar pho-tosphere at R ∗ = r ( τ = 2 / , we provide unique mass-lossrates.As shown in Sect. 2 this method provides a fully an-alytic mass-loss rate for an isothermal wind. Using a gridof numerical models, we further provide a fit to the cor-rection factor between the isothermal and non-isothermalmass-loss rate for the same input parameters (Sect. 3.2).Combining this correction factor with the isothermal mass-loss rates then presents a fully analytic mass-loss rate pre-scription for red supergiants, given by the combination ofEquations 5, 13, 14, and 25. Comparisons of this mass-lossrate with both observations and empirical mass-loss rateprescriptions show that such a turbulent pressure drivenmass-loss rate is able to recover realistic mass-loss ratesover the range of parameters appropriate to red super-giants. These comparisons also inform our suggested choiceof v turb = 18 . ± . km s − as an appropriate representa-tive value of v turb in mass-loss rate calculations where thisparameter is not well constrained.To continue the discussion, we enumerate some itemsof future work required to place this analytic, theoreticalmodel of red supergiant mass-loss on even firmer footing.From an observational side, as an immediate first item es-timations of turbulent velocities for a wider body of redsupergiants are needed. As it stands, the limited samplethat we have examined here provides highly encouragingresults, and additional objects would allow the model to becalibrated more carefully. Additionally, such observations would be of key importance to understand how this currentfree parameter of the model scales with stellar parameters.Finally, such future observations could also help constrainthe depth dependence of v turb . Prior work has suggestedthat turbulent velocity may actually increase with distancefrom the star over at least part of the wind launching re-gion (de Koter et al. 1988; Josselin & Plez 2007), and betterunderstanding this profile would be highly useful in furtherdeveloping the model presented here.To further develop the model itself, main efforts mayfocus on relaxing key approximations made in this paper.One of the first of these would be to use a radial profileof v turb , perhaps as inferred by observations above, ratherthan the spatial average employed here. Such a radial profilemight arise from dissipation of turbulence by shocks. Inclu-sion of a radial profile in v turb would allow, for instance,for a prediction of the terminal velocities of red supergiantwinds instead of only mass-loss rates, as the energy avail-able to drive the wind would no longer diverge with in-creasing outer radius of the model. Further generalizationof v turb could also include variations over the surface of thestar and/or in time. Both of these would arise naturally ifturbulence is related to surface convection in the star, andtheir inclusion would help to explain the observed clumpystructure of the wind.A further investigation of the overall structure of red su-pergiant winds may also re-introduce dust and/or moleculesand their associated opacity to examine the impact thiswould have for the proposed model. As discussed in the in-troduction, previous models have assumed that radiationpressure is the dominant driving force in carrying the massloss of red supergiants, with the main question being theorigin of the “missing force” necessary to get material fromthe stellar photosphere up to the region where moleculesand dust can form, the latter occurring potentially sev-eral stellar radii away from the photosphere. Seemingly cir-cumventing such a model, we show that turbulent pressureextended red supergiant atmospheres are already able togenerate the appropriate mass-loss rates in the absence ofradiation pressure on dust or molecules. However, it is im-portant to recognize that this does not necessarily meanthat enhanced opacity from dust or molecules has a negli-gible contribution. In fact, these opacities almost certainlyplay a crucial role in setting the wind terminal speed as tur-bulent pressure dies off away from the stellar surface (seethe discussion in Sect. 3.1). It is also clear, as discussedin Sect. 2.2, that the overall scale of opacity and its radialprofile may contribute to the optical depth and scale heightof the wind and thus to the mass loss rate retrieved by thismodel. Further, it is even possible that the radial profile ofturbulent pressure is such that turbulent pressure plays thedominant role in levitating material off the stellar surface,but that dust or molecular opacity carries the wind throughthe critical point as previously theorized.Another important line of investigation regards the ori-gin of the turbulent velocities and extended atmospheresthis model leverages. As discussed by Arroyo-Torres et al.(2015), previous models have focused on convection through3D, so-called “star-in-a-box” radiation hydrodynamics sim-ulations such as those computed with CO BOLD (Frey-tag et al. 2012), and on pulsations through 1D simulationsakin to those used for AGB stars and especially Mira vari-abls (e.g.
CODEX
Ireland et al. 2008, 2011). While convec-tion and pulsations in the stellar atmosphere both provide
Article number, page 11 of 14 &A proofs: manuscript no. rsg_modParker_final tempting options, convection simulations generally fail togenerate the turbulent velocities or atmospheric extensionsinferred from observations, and pulsations are required tobe of unrealistically high amplitude compared to observedlightcurves for red supergiants (see, e.g. Josselin & Plez2007; Arroyo-Torres et al. 2015).Similarly, one can examine the role that rotation mighthave in generating similar levels of extension to the stellaratmosphere. This can be done by applying the centrifu-gal term ( Ω r ) either along side or instead of v turb in thepreceding analysis. Taking the limiting case where turbu-lent pressure is omitted, we find that matching the atmo-spheric extension of a 18.2 kms − turbulent velocity even atthe most optimal position on the stellar equator (includingthe effects of stellar oblateness) requires a rotation veloc-ity almost of the orbital velocity, implausibly fast fora red supergiant. As we reinforce through this paper thatatmospheric extension is crucial for the mass-loss process,further examination of the possible missing ingredients forsuch models as well as potential constructive interactionsbetween them is now highly timely.To wrap up the discussion, we address some of the po-tential implications of applying our mass-loss rates, espe-cially for the case of stellar evolution modeling. As discussedin Sect. 4.2, the current state-of-the-art mass-loss rates forstellar evolution modeling are all empirical. Given that noneof the empirical rates match the theoretical one we presenthere for all stars, these new theoretical rates should impactthe mass, luminosity, and radius distributions of supergiantsupernova progenitors. Additionally, the steeper trend ofincreasing mass-loss rate with increasing luminosity wouldnaturally turn more massive stars evolving toward the redsupergiant branch back to the blue due to envelope strip-ping, thereby predicting a decreased upper luminosity limitfor red supergiants. Finally, these rates provide intriguingimplications for stellar evolution in lower metallicity en-vironments. Many standard stellar evolution implementa-tions of mass loss in the red supergiant branch impose adownturn with metallicity partially inspired by the ˙ M ∝ Z n dependence of line-driven hot star winds (Vink et al. 2001;Mokiem et al. 2007; Björklund et al. 2020). Observationsanalyzed by Mauron & Josselin (2011) do seem to suggestsuch a scaling of ˙ M ∝ ( Z/Z (cid:12) ) . from the Galaxy to the redsupergiants in the Small Magellanic Cloud. However, theyconcluded that the scaling to the Large Magellanic Cloudwas not well constrained by the samples they considered,emphasizing that more investigation is needed. Meanwhile,the model we present here should not be strongly dependenton metallicity. While the optical depth and scale height ofthe wind clearly do depend on dust and molecular content,which in turn depend on metallicity, this is a weaker sec-ondary scaling when compared with the turbulent pressureitself. Moreover, decreasing opacity can lead to increasedmass-loss rates in this model depending where in parame-ter space a star sits (see Sect. 2.2). So long as the turbulentpressure mechanism itself is not strongly metallicity depen-dant, which turbulence triggered by hydrogen recombina-tion is not, this could then imply much stronger mass losson the red supergiant branch at low metallicity than pre-viously assumed, and as such significantly different stellarevolution tracks. Acknowledgements.
We would like to thank the anonymous refereefor their comments which enhanced the discussion in this paper. We acknowledge support from the KU Leuven C1 grant MAESTROC16/17/007. L.D. also acknowledges funding from the European Re-search Council (ERC) under the European Union’s Horizon 2020research and innovation programme (grant agreements No. 646758:AEROSOL with PI L. Decin).
References
Arroyo-Torres, B., Wittkowski, M., Chiavassa, A., et al. 2015, A&A,575, A50Barber, R. J., Tennyson, J., Harris, G. J., & Tolchenov, R. N. 2006,MNRAS, 368, 1087Beasor, E. R., Davies, B., Smith, N., et al. 2020, MNRAS, 492, 5994Björklund, R., Sundqvist, J. O., Puls, J., & Najarro, F. 2020, arXive-prints, arXiv:2008.06066Bladh, S. & Höfner, S. 2012, A&A, 546, A76Bladh, S., Höfner, S., Nowotny, W., Aringer, B., & Eriksson, K. 2013,A&A, 553, A20Braun, K., Baade, R., Reimers, D., & Hagen, H. J. 2012, A&A, 546,A3Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115De Beck, E., Decin, L., de Koter, A., et al. 2010, A&A, 523, A18de Jager, C., Nieuwenhuijzen, H., & van der Hucht, K. A. 1988, A&AS,72, 259de Koter, A., de Jager, C., & Nieuwenhuijzen, H. 1988, A&A, 200,146Decin, L., Hony, S., de Koter, A., et al. 2006, A&A, 456, 549Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146Freytag, B., Steffen, M., Ludwig, H. G., et al. 2012, Journal of Com-putational Physics, 231, 919Friend, D. B. & Abbott, D. C. 1986, ApJ, 311, 701Goorvitch, D. 1994, ApJS, 95, 535Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486,951Gustafsson, B. & Plez, B. 1992, in Instabilities in Evolved Super- andHypergiants, ed. C. de Jager & H. Nieuwenhuijzen, 86Harwit, M., Malfait, K., Decin, L., et al. 2001, ApJ, 557, 844Hearn, A. G. 1975, A&A, 40, 355Heger, A. & Langer, N. 1998, A&A, 334, 210Heger, A. & Langer, N. 2000, ApJ, 544, 1016Höfner, S. & Olofsson, H. 2018, A&A Rev., 26, 1Iglesias, C. A. & Rogers, F. J. 1996, ApJ, 464, 943Ireland, M. J., Scholz, M., & Wood, P. R. 2008, MNRAS, 391, 1994Ireland, M. J., Scholz, M., & Wood, P. R. 2011, MNRAS, 418, 114Josselin, E. & Plez, B. 2007, A&A, 469, 671Langer, N. 1997, in Astronomical Society of the Pacific ConferenceSeries, Vol. 120, Luminous Blue Variables: Massive Stars in Tran-sition, ed. A. Nota & H. Lamers, 83Levesque, E. M. 2017, Astrophysics of Red Supergiants (IOP Publish-ing)Lucy, L. B. 1971, The Astrophysical Journal, 163, 95Maeder, A. & Meynet, G. 2000, A&A, 361, 159Mauron, N. & Josselin, E. 2011, A&A, 526, A156McKemmish, L. K., Masseron, T., Hoeijmakers, H. J., et al. 2019,MNRAS, 488, 2836–2854Meynet, G. & Maeder, A. 2000, A&A, 361, 101Mokiem, M. R., de Koter, A., Vink, J. S., et al. 2007, A&A, 473, 603Nieuwenhuijzen, H. & de Jager, C. 1990, A&A, 231, 134O’Gorman, E., Harper, G. M., Ohnaka, K., et al. 2020, A&A, 638,A65Ohnaka, K., Hofmann, K. H., Schertl, D., et al. 2013, A&A, 555, A24Ohnaka, K., Weigelt, G., & Hofmann, K. H. 2017, Nature, 548, 310Owocki, S. P., Cranmer, S. R., & Gayley, K. G. 1996, ApJ, 472, L115Parker, E. N. 1958, ApJ, 128, 664Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10Petrenz, P. & Puls, J. 2000, A&A, 358, 956Reimers, D. 1975, Memoires of the Societe Royale des Sciences deLiege, 8, 369Reimers, D. 1977, A&A, 61, 217Rogerson, J. B., J. & Lamers, H. J. G. L. M. 1975, Nature, 256, 190Schrijver, C. J. & Zwaan, C. 2000, Solar and Stellar Magnetic Activity(Cambridge University Press)Sobolev, V. V. 1960, Moving envelopes of stars (Harvard UniversityPress)Sylvester, R. J., Skinner, C. J., & Barlow, M. J. 1998, MNRAS, 301,1083van Loon, J. T., Cioni, M. R. L., Zijlstra, A. A., & Loup, C. 2005,A&A, 438, 273van Loon, J. T., Groenewegen, M. A. T., de Koter, A., et al. 1999,A&A, 351, 559Verhoelst, T., van der Zypen, N., Hony, S., et al. 2009, A&A, 498, 127Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369,574Wood, P. R., Bessell, M. S., & Fox, M. W. 1983, ApJ, 272, 99
Article number, page 12 of 14. D. Kee et al.: Analytic RSG mass-loss rates from turbulent pressure
Appendix A: Computation of g rad and κ fromaccelerating molecules As discussed in Sect. 2, the total opacity scale of the windswe treat here can be strongly impacted by the Doppler shift-ing of spectral lines, leading to opacities exceeding the staticRosseland mean opacity. To show this, we begin with a gen-eralized form of radiative acceleration in the ˆ n direction, g rad = 1 c (cid:73) (cid:90) ∞ κ ν I ν dν ˆn d Ω , (A.1)with speed of light c , extinction coefficient κ ν [ cm g − ],and intensity per unit frequency I ν . If one assumes that theextinction comes from a single, isolated spectral line, thena natural choice is to split κ ν into a normalized, frequency-dependent shape of the spectral line, or profile function φ ν ,and a line-integrated total extinction coefficient κ L = σ cl n l f lu ρ (cid:18) − exp (cid:18) − h ν o k B T (cid:19)(cid:19) . (A.2)In this expression for κ L , σ cl is the classical oscillator crosssection, n l is the number density of particles in the lowerlevel of the transition, f lu is the quantum mechanical os-cillator strength of the transition, and ρ is mass density.The final exponential term accounts for stimulated emissionunder the assumption of local thermodynamic equilibrium(LTE). Such stimulated emission can play an importantrole for the considered infrared spectral lines transitions asthe energy of the transition, given by Planck’s constant h times rest frequency ν o , can be comparable to the thermalenergy of the gas, given by Boltzmann’s constant k B timesgas temperature T . Note that here κ L has picked up unitsof [ cm g − Hz ] as the profile function is defined per unitfrequency. Under these substitutions, the equation we nowneed to solve becomes g line = κ L c (cid:73) (cid:90) ∞ φ ν I ν dν ˆn d Ω . (A.3)At this point, it is important to recognize that the in-tensity fed into this expression can itself be a strong func-tion of frequency even in the presence of an optically thincontinuum, as is the case for red supergiants, due to theattenuation of intensity by the spectral line we consider it-self. Therefore, we have to obtain intensity from the formalsolution of the radiative transfer equation I ν = I o ν e − τ ν + (cid:90) S ν ( t ν ) e −| τ ν − t ν | dt ν . (A.4)While in general this becomes quite complex as the sourcefunction S ν depends on I ν , we can simplify by leveragingsymmetry arguments. In the Sobolev limit considered here(see further below), the second term here becomes fore-aftsymmetric such that it cancels out in the integral (cid:72) ˆ nd Ω .As regards the optical depth τ ν , we take the approximationthat this arises only from extinction from the spectral lineitself such that τ ν = (cid:90) κ L φ ν ρ dl , (A.5) where the mixed frequency and spatial notations can berectified by using the Doppler formula ν − ν o ν o = v l c , (A.6)such that τ ν = (cid:90) κ L φ ν ρ (cid:18) ∂l∂v l (cid:19) (cid:18) ∂v l ∂ν (cid:19) dν = (cid:90) κ L ρ cν o ∂v l /∂l φ ν dν . (A.7)Finally, we can take the Sobolev approximation (Sobolev1960) which argues that all hydrodynamic variables are ef-fectively constant over a line resonance region and thereforecan be pulled through the frequency integral to give g line = κ L c (cid:73) (cid:90) ∞ φ ν I o ν e − τ S (cid:82) ∞ φ ν dν dν ˆn d Ω , (A.8)where we have introduced the Sobolev optical depth τ S = κ L ρc/ ( ν o ∂v l /∂l ) . This final expression can now beanalytically solved to give g line = κ L c (cid:73) I o ν (cid:18) − e − τ S τ S (cid:19) d Ω , (A.9)which further simplifies in the case of a spherically sym-metric flow and radiation from a point star to g line , r = π κ L I o ν c (cid:18) − e − τ S,r τ S,r (cid:19) , (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 employedhere is valid in wind outflow regions where the flow speedexceeds a few times the characteristic thermal speed, v th ,such that there is not significant absorption of intensity en-tering the resonance region by a spectral line in the hy-drostatic stellar atmosphere, and such that the physicalsize of the resonance region is small compared to the scaleover which the hydrodynamic quantities entering τ S vary.In practice, this condition is nearly always met, however, asthe thermal speeds of the molecules we examine are reducedcompared to the bulk wind sound speed by the potentiallyquite substantial ratio of mean molecular mass of the bulkgas to mass of the molecule. Therefore, as this expressionholds for each individual line over the bulk of the wind,to generalize this to a spectrum of lines one simply needsto sum g line . Numerous tabulated line lists are available tofacilitate this process, each including ν o , either f lu or theEinstein coefficients necessary to compute it, and the energyand degeneracy of the levels involved in each of the transi-tions, which can be used to compute n l , for instance fromthe Boltzmann distribution in Local Thermodynamic Equi-librium. Taking this fully summed acceleration, it is thenstraightforward to back out a flux-weighted mean opacity κ of the line acceleration by appealing to Equation 3 suchthat κ = cF r (cid:88) g line , r , (A.11) Article number, page 13 of 14 &A proofs: manuscript no. rsg_modParker_final where the summation is over all lines. As more force is avail-able from an accelerating spectral line than a static one,this will always return flux-weighted mean opacity that islarger than the Rosseland mean opacity, thereby substanti-ating our choice to test the effects of varying κ in Sect. 2.2.in Sect. 2.2.