Boil-off of red supergiants: mass loss and type II-P supernovae

The mass loss mechanism of red supergiant stars is not well understood, even though it has crucial consequences for their stellar evolution and the appearance of supernovae that occur upon core-collapse. We argue that outgoing shock waves launched near the photosphere can support a dense chromosphere between the star's surface and the dust formation radius at several stellar radii. We derive analytic expressions for the time-averaged density profile of the chromosphere, and we use these to estimate mass loss rates due to winds launched by radiation pressure at the dust formation radius. These mass loss rates are similar to recent observations, possibly explaining the upward kink in mass loss rates of luminous red supergiants. Our models predict that low-mass red supergiants lose less mass than commonly assumed, while high-mass red supergiants lose more. The chromospheric mass of our models is $\sim$0.01 solar masses, most of which lies within a few stellar radii. This can help explain the early light curves and spectra of type-II P supernovae without requiring extreme pre-supernova mass loss. We discuss implications for stellar evolution, type II-P supernovae, SN 2023ixf, and Betelgeuse.

1. INTRODUCTION It has long been known that red giants and red supergiants (RSGs) lose mass through dense, slow-moving winds (e.g., Reimers 1975).These winds can carry mass away at high rates 10 −8 M ⊙ /yr ≲ Ṁ ≲ 10 −4 M ⊙ /yr (e.g., Antoniadis et al. 2024), partially or fully stripping the hydrogen envelopes of RSGs before their ultimate death via a supernova explosion or collapse to a black hole.RSGs are formed from massive stars (10 M ⊙ ≲ M ≲ 40 M ⊙ ) during their post-main sequence evolution, with large radii of 300 R ⊙ ≲ R ≲ 1000 R ⊙ , high luminosities of 10 4 L ⊙ ≲ L ≲ 3 × 10 5 L ⊙ , and cool surface temperatures of 3500 K ≲ T eff ≲ 4500 K (e.g., Levesque et al. 2005).Their low escape speeds v esc ∼ 100 km/s and high luminosities probably contribute to their high mass loss rates, yet the physical mechanisms underlying mass loss are not well understood.
It is widely accepted that these winds are accelerated by radiation pressure on dust (e.g., Hoyle and Wickramasinghe 1962), which forms well above their photospheres where the temperature drops to T ≲ 1500 K (Field 1974;Höfner and Olofsson 2018) such that dust can form (Figure 1).The high opacity of dust means that dust-laden gas has a large enough opacity to exceed the Eddington limit, such that it will be accelerated outwards by the momentum of photons emitted from the star.However, the dust can only form at radii several times larger than the photosphere of the star, R d ≳ 5 R, jfuller@caltech.eduand the mechanisms by which material is lofted up to the dust formation radius are highly uncertain.This windlaunching problem for RSGs is crucial for determining their mass loss rates, and it is the main focus of this paper.
Perhaps the most detailed investigations of mass loss from cool stars have been performed for asymptotic giant branch (AGB) stars (see Höfner and Olofsson 2018 for a review).Although these stars have lower mass and smaller luminosities than RSGs, they have similar surface temperatures and escape speeds, and the mass loss mechanism is believed to be similar to that described above.For AGB stars, some combination of turbulent pressure and shock waves driven by radial pulsations are thought to propel gas above the photosphere (Höfner et al. 2003).Because this chromospheric material is dense, it can cool (unlike the Sun's corona, Cranmer and Saar 2011) so that its temperature maintains radiative equilibrium with the star and decreases approximately as T ∝ r −1/2 in the chromosphere.Evermore sophisticated radiation hydrodynamics simulations of this process (Freytag and Höfner 2008;Freytag et al. 2017;Freytag and Höfner 2023) support this general picture and also give detailed predictions for atomic, molecular, and dust spectral features in the spectra of these stars.
For the last decade, transient surveys have been detecting type II-P supernovae at increasingly early times (sometimes just hours after shock-breakout) when the appearance of the supernova is very sensitive to the properties of the exploding RSG and its immediate sur-roundings (e.g., Yaron et al. 2017;Hiramatsu et al. 2021;Bruch et al. 2021;Li et al. 2023;Meza Retamal et al. 2024).With typical ejecta velocities of ∼ 10 9 cm/s, the outgoing SN ejecta takes ∼5 days to expand from the star's photosphere out to the dust formation radius of R d ∼ 10 R ∼ 5 × 10 14 cm.Early SN observations thus probe the SN as it is sweeping up material between the star's surface and dust formation radius.The results of these observations have been clear: a large fraction of type II-P SNe show indications of interaction between the SN ejecta and dense circumstellar material (CSM) within the dust formation radius (Khazov et al. 2016;Bruch et al. 2021Bruch et al. , 2023)).
This CSM is orders of magnitude denser than expected from stellar evolution models of RSGs, which predict sharp drops in density above the photosphere.The CSM is also a few orders of magnitude denser than expected if the CSM is formed by a wind with mass loss rate and velocity typical for RSGs, Ṁ ∼ 10 −6 M ⊙ /yr, v wind ∼ 30 km/s, (Dessart et al. 2017;Moriya et al. 2017;Morozova et al. 2017;Moriya et al. 2018;Morozova et al. 2018).The frequent conclusion of SN observers and modelers has thus been that RSGs lose mass at much higher rates during the final years of their lives.This idea has been supported by some RSG models that have predicted pre-supernova outbursts that eject mass in the weeksmonths before explosion (Woosley and Heger 2015;Fuller 2017).
In this paper, we argue that the RSG wind-launching problem and the type-II SN interaction problem are inextricably linked.We show that CSM surrounding exploding RSGs is naturally produced by a dense, radially extended chromosphere between the stellar surface and dust formation radius.This chromosphere is approximately in hydrostatic equilibrium in a time-averaged sense, and it is supported by the momentum flux of upward-traveling weak shock waves (analogous to turbulent pressure support) excited by vigorous convection near the stellar photosphere.By computing the density of the chromospheric material at the dust formation radius, we also predict RSG mass loss rates.We show that these models can explain several features of RSG mass loss rates as a function of RSG luminosity.Finally, we show that the extended chromospheres of these RSGs helps to explain the unexpected early light curves and flash-ionized emission lines observed in many type II-P SNe.

CHROMOSPHERIC DYNAMICS
2.1.Equations of Motion In red supergiants, both 1D stellar models with convective mixing length theory (MLT) and 3D stellar models predict typical convective speeds near the photosphere of v con ∼ 10 km/s (Goldberg et al. 2022a).In comparison, the escape speeds from the are where M and R are the mass and photospheric radius of the star.Hence, typical convective elements move slower than the escape speed and only rise a small height above

Dust Formation Zone
Wind Acceleration Shock-supported Chromosphere Fig. 1.-Cartoon illustrating the chromospheric and wind dynamics of red supergiants (not to scale).Outgoing shock waves launched by vigorous convective motions support optically thin material above the photosphere, forming an extended and dynamically varying chromosphere.Material that rises above ∼5 stellar radii forms dust, increasing its opacity and creating a radiationdriven wind.
Simulations show that turbulent pressure support dominates over thermal pressure support for material above the photosphere.For a turbulent pressure P turb ∼ ρv 2 turb , where ρ is the density and v turb is the RMS turbulent velocity, we naively expect a chromospheric scale height of H ∼ P turb /(ρg) ∼ v 2 con /g ∼ 2r(v turb /v esc ) 2 .This is a factor of ∼50 smaller than the stellar radius for typical RSGs, and so we might expect the density to fall very steeply if it is supported by turbulent pressure.
However, an important factor to consider is the intermittency of convection, which causes large fluctuations in the convective velocities near the surface of the star.Some regions of the star, or some epochs in time, have larger convective velocities and will have larger chromospheric scale heights.We show this can greatly increase the chromospheric density and associated mass loss rate.
To build a model for the chromospheric structure, we consider hydrostatic equilibrium between gravity and momentum depostion from outgoing shock waves.The momentum equation can be written in the form The ρ⃗ v term is proportional to the mass flux, which may be non-zero but is assumed to be constant in space/time.Hence, averaging over time, this term vanishes.For simplicity, we will neglect thermal pressure since it is small for moderate or strong shocks (see Bertschinger and Chevalier 1985) for a model including thermal pressure).We also neglect radiative forces because the low opacity and density of the chromosphere causes it to be optically thin, but later we will account for radiative acceleration of dust grains in the dust formation region.Rather than assuming isotropic turbulent motion, we consider purely radial motion, which is a better model to capture the behavior of shocks propagating nearly radially out into the chromosphere.Shocks in RSG chromospheres are nearly isothermal due to the short radiative relaxation time in the post-shock gas (Schirrmacher et al. 2003).Hence, to good approximation we expect an adiabatic index of γ ≃ 1.For a strong nearly isothermal shock, the post-shock and post-cooling velocity v of a fluid element is nearly equal to the shock velocity v s (e.g., Shu 1992).Hence, the radial momentum flux of an outward going shock wave is ρv 2 s , where v s is the shock speed.
Applying these simplifications to equation 2, hydrostatic equilbrium requires where r is the radial coordinate.Here, ρ should be regarded as the time-averaged density, whereas the actual density and fluid velocity fluctuate in time as shock waves pass by.We approximate M as constant above the photosphere because the chromospheric mass is a very small fraction of the star's total mass.Simulations (e.g., Freytag et al. 2017;Goldberg et al. 2022a) show that the turbulent velocity is nearly constant in the region above the photosphere of a AGB or RSG star.This may be explained by the steepening of turbulent motions into moderate shocks (Mach numbers M ∼ 3) in the chromosphere.Such moderate shocks propagating down a density gradient typically steepen until their Mach numbers are of order unity, and then maintain nearly constant Mach number (Coughlin et al. 2018) and hence nearly constant speed in the nearly isothermal chromosphere.
Another justification for a nearly constant shock speed above the photosphere is the isothermal nature of the shocks produced.Following Matzner and McKee (1999), a shock wave propagating through a star accelerates the material to a velocity of Here, v 0 , m 0 , ρ 0 , and r 0 are the shock velocity, mass coordinate, density, and radius at an arbitrary reference location.The value of β depends on the equation of state of the gas (Sakurai 1960), but is approximately β ≈ [2 + 2γ/(γ − 1)] −1 for strong shocks (Whitham 2011), which gives µ ≃ 0.2 for γ = 4/3 as appropriate for adiabatic and radiation-pressure dominated post-shock gas within stellar interiors.But for nearly isothermal shocks appropriate for RSG chromospheres, γ ≃ 1 and β ≃ 0. The mass coordinate m is nearly constant for a chromosphere with small total mass.In this case, from equation 4, the post-shock velocity v s ∼ v 0 is nearly constant as the shock propagates outward.
A third justification for a constant shock speed as a function of radius are the models of Bertschinger and Chevalier (1985), who examined atmospheres of Mira variables supported by a periodic train of outgoing shocks.In the plane-parallel limit, they showed that periodic solutions have shock velocities that are constant as a function of height, while the density falls exponentially (when averaging over time), as shown by the calculation below.The plane-parallel limit applies when the separation between subsequent shocks is less than the stellar radius, which is only marginally true in our case.Simulations (e.g., Goldberg et al. 2022a) show that shocks are launched approximately on the convective turnover time t con ∼ 0.5R/v con such that the separation between shocks is ∼ 0.5R.
The shock support term can also be derived heuristically.For strong shocks passing through a medium, the rate of momentum deposition per unit radius is ṗ = d(4πr 2 ρv s )/dt.For a time-steady state, the time derivative is d/dt → v s d/dr.If v s is constant as a function of radius, it can be moved inside the derivative.The momentum deposition per unit volume is thus ṗ/(4πr 2 ) = (1/r 2 )d(ρr 2 v 2 s )/dr.Remarkably, in the steady-state model, the momentum deposition rate is independent of the frequency at which shocks pass by a given point (i.e., number of shocks per unit time), and is only dependent on the shock amplitude v s .

Chromospheric Structure
Proceeding with the assumption that v s is constant as a function of radius, equation 3 can be solved to yield the density profile This is similar to the structure of an isothermal atmosphere, but with the sound speed replaced by v s , and with the extra geometric factor of (R/r) 2 due to the use of the divergence of the momentum flux rather than a turbulent pressure.Just above the stellar surface, equation 5 is identical to the density profile found by Bertschinger and Chevalier (1985) for a periodic train of shocks in a plane-parallel atmosphere.Importantly, the density profile depends only on the amplitude of shocks, and not on their temporal frequency.
As mentioned above, shocks are generated by turbulent convection so that v s is not actually constant, but fluctuates stochastically in time.Since this term enters into the exponential of equation 5, this could cause huge changes in the chromospheric density profile.We assume the convective velocity distribution follows a Gaussian probability distribution where v con is the typical convective velocity.From this definition, the mean turbulent velocity is v con .We use a Gaussian rather than Maxwellian distrubtion because we consider only the radial component of the convective velocities.
If we assume that v s is constant throughout the chromosphere, the time-averaged density profile is 1−R/r vcon .
(7) Remarkably, accounting for fluctuating convective velocities transforms the isothermal atmosphere density profile of equation 5 into an exponential distribution.This is important because the ratio of v esc /v con ∼ 10 is a large number, so the time-averaged density profile of equation 7 falls off much more slowly than that of equation 5.The integrand of equation 7 is maximized where v s ∼ √ v esc v con ∼ 3v con , hence the time-averaged density profile at large radii is dominated by uncommon moments of more intense turbulent motion where the convective velocities are ∼2 sigma larger than their usual values.
Figure 2 shows the density profiles of equation 5 and equation 7 above a RSG photosphere.The underlying model is a 13.8 M ⊙ RSG with R = 844R ⊙ generated using MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015)).This model had a ZAMS mass of 15 M ⊙ and has been evolved to core-collapse.We use MLT with α MLT = 2.5 to measure the convective velocity, which we evaluate where the optical depth is τ = P rad c/(P v con ), where P rad and P are the radiation and total pressures.This yields approximately correct RMS convective velocities (see Goldberg et al. 2022a), which for RSGs are quite similar to the Time-averaged -Top: Profile of overlying chromospheric mass as a function of radius corresponding to the models in Figure 2. We expect chromospheres of ∼ 10 −2 − 10 −1 M ⊙ , mostly within ∼2 stellar radii.Bottom: Time-averaged mass profiles of the same model.sound speed c s .This model has v con ≃ 0.08v esc near the star's surface.
We plot the density profiles of equation 5 for v s /v esc =0.08, 0.16, and 0.24 to demonstrate how moments of large v s greatly increase the chromospheric density.We also plot the time-averaged density profile of equation 7 for values of v con /v esc =0.06, 0.08, and 0.10.These values are representative of our RSG models, with low-mass RSGS at the lower end of this range and highmass RSGS or pre-SN RSGs at the upper end.For this 15 M ⊙ model, the corresponding range in mass loss rate is Ṁ ∼ 10 −7 − 10 −4 M ⊙ /yr (see Section 3), demonstrating its senstivity to v con /v esc .
Figure 2 shows that the chromospheric density falls off steeply within one stellar radius above the photosphere, but then approaches a r −2 profile.Nonetheless, within a few stellar radii, the chromospheric density is typically orders of magnitude larger than that implied by simple constant-velocity wind models.For comparison, a constant-velocity wind at mass loss rate Ṁwind has a density We plot ρ wind for v wind = 30 km/s and three different mass loss rates.The chromospheric density profiles of our models can be much larger than those of simple wind models below a few stellar radii, even for large mass loss rates of Ṁwind = 10 −4 M ⊙ /yr.Between 1-2 stellar radii, our chromospheres may have densities comparable to an extreme wind of Ṁwind = 10 −2 M ⊙ /yr, even though the true mass loss rates of our models are orders of magnitude smaller (Section 3).
Correspondingly, the circumstellar mass is much larger in our models than naive expectations.We define the overlying CSM mass as where R d is the dust formation radius.Figure 3 shows the overlying CSM mass for both instantaneous and time-averaged cases.We find time-averaged chromospheric masses out to R d (see next section) of ∼ 10 −2 M ⊙ , with large variation for different values of v s /v esc or v con /v esc .Most of this mass is near the photosphere, but the time-averaged mass above two stellar radii can be ∼ 10 −3 M ⊙ .For comparison, a constant wind model with a similar mass loss rate of Ṁ = 10 −5 M ⊙ /yr contains only ∼ 10 −5 M ⊙ within two stellar radii.

MASS LOSS 3.1. Dust Formation and Wind Acceleration
If chromospheric material extends to large radii, it can cool enough to form large amounts of dust.Radiation -The velocity profiles of the chromosphere/wind for the same models as Figure 4, showing how the wind speeds are small below the dust formation radius, but approach the terminal wind speed at large distances.
pressure on this dust will then accelerate this material outwards, driving mass loss from the RSG.Except just below shocks, the temperature profile of the chromosphere is approximately in radiative equilibrium with the star due to molecular and dust absorption/emission processes (Schirrmacher et al. 2003).The grey radiative equilibrium temperature for an optically thin chromosphere in the Eddington approximation is (e.g., Chandrasekhar 1934;Bowen 1988;Winters et al. 1997) where T * is the photospheric temperature.For a given dust formation temperature T d , the corresponding dust formation radius is For RSG surface temperatures of T * ≈ 3600 K and dust formation temperatures of T d ∼ 1200 K, dust condenses at a dust formation radius of R d ∼ 4.5 R.
To estimate the wind velocity profile due to radiative driving of dust in our model, we approximate the radiative acceleration as zero below R d and falling off as 1/r 2 above the dust formation radius.We also neglect turbulent pressure support above R d so that the momentum equation in this region is Here, the constant f Ed = L/L Ed determines the magnitude of the radiation force on the dust-gas mixture and will be computed below.For an initial velocity v(R d ) = v d at the dust formation radius, the solution to this equation is where To match our solution above and below the dust formation radius, we choose an outflow velocity v d = v con at the dust formation radius.This is motivated by the fact that the dust formation radius is where the flow transitions from an approximately hydrostatic to an unbound outflow, such that the outflow velocity begins to exceed v con , and the advection term becomes important in the momentum equation.The mass loss rate is then This procedure is similar to computing the mass loss of a Parker wind at the sonic point.However, in our problem the sonic point is not important because the dust formation radius is at much smaller radius than the sonic radius (i.e., the effective Bondi radius), which lies at R s = GM/v 2 con ∼ 50R.The mass loss rate is determined primarily by the chromospheric density at the dust formation radius, and it is only weakly sensitive to the radiative acceleration of the wind above this point.
We choose to use observed wind velocities v ∞ to constrain x d , but one could use dust opacities and mass fractions to calculate x d and then predict the corresponding value of v ∞ .From equation 14, we have In this work, we use v ∞ = 30 km/s, in line with observational estimates (Sjouwerman et al. 1998;van Loon et al. 2001;Marshall et al. 2004), although this velocity likely depends weakly on the stellar luminosity (e.g., Mauron and Josselin 2011).This requires x d ∼ 2, approximately what we expect for realistic dust-laden gas opacities of κ d ∼ 5 cm 2 /g (Höfner et al. 2003).
To compute the time-averaged outflow speed below the dust formation radius, we simply require Ṁ to be constant with radius such that v out = Ṁ /(4πρr 2 ).Finally, to compute the time-average density above the dust formation radius, we again use the constant mass loss rate so that the density is (16) As mentioned above, dust formation begins at temperatures ranging from T ≃ 1000 − 2000 K, with different species condensing at different temperatures.Here we set the dust formation radius to R d = 5R, a reasonable value for RSGs.Different choices for R d would result in slightly different density/velocity profiles above R d .The mass loss rate of our model is only weakly sensitive to the value of R d because the factor of R −2 d in the density profile (equation 5) cancels out the factor of R 2 d in the mass loss rate (equation 15).In reality, dust starts to condense over a range of temperatures/radii, and there is no single dust formation radius.
Figures 4 and 5 shows our computed time-averaged density and velocity profiles of the star's CSM for different values of v con /v esc .The density profiles are the same as Figure 2 below the dust formation radius, but Fig. 6.-Mass loss rates of stellar models as a function of their luminosity.The small gray symbols are measured mass loss rates of red supergiants in the LMC from Antoniadis et al. 2024, while colored symbols are our models.We only show the red supergiant phase of evolution, when T eff < 4200 K.We plot one symbol for every 10 4 years of evolution, and use a large star symbol to denote the pre-supernova mass loss rate.The full evolution of 30 M ⊙ and 40 M ⊙ models is not shown because they turn into yellow supergiants after heavy mass loss.An alternative mass loss mechanism is required to explain stars with L ≲ 2 × 10 4 L ⊙ .
we now extend the plot above R d .The velocity profile has low outflow velocities below R d , before the transition into a wind above the dust formation radius.The density is much larger than that of a constant-velocity wind below R d , and then it approaches the constant-velocity solution above R d where the wind speed approaches is asymptotic value.
For the time-averaged models, there is a downward kink in the density profiles at the dust formation radius R d .This is caused by the upward kink in the outflow velocity at R d (Figure 5) due to radiative acceleration.While the true velocity and density profiles are unlikely to contain sharp kinks, we do expect a small dip in CSM density above dust formation radius due to the acceleration of the wind.

Red Supergiant Mass Loss Rates
We implement the mass loss rate of the previous section (equation 15) into stellar evolution models generated with MESA as described above.At surface temperatures above 4500 K we use the "Dutch" wind prescription with scaling factor η = 0.5.At cooler surface temperatures, we use our wind prescription.Our MESA inlists and run_star_extras files are available at this link.
Figure 6 shows the predicted mass loss rates for RSGs with masses of 10 − 40 M ⊙ .We only plot the RSG phase where T eff < 4200 K, where the surface convection zone is deep.Usually the mass loss rate from boil-off is negligible at higher temperatures when the star is more compact.However, those models can have narrow super-Eddington sub-surface convective zones and we do not trust the MLT prediction for v con in those models, so we do not investigate those cases here.
Our mass loss rates are similar for a given luminosity to those recently measured for RSGs in the LMC (Antoniadis et al. 2024;Wen et al. 2024).These mass loss rates are substantially smaller than those measured in the SMC (Yang et al. 2023), but those may be systematically overestimated according to Antoniadis et al. 2024.Most importantly, the predicted scaling of mass loss with luminosity above L > 3 × 10 4 L ⊙ is quite similar between our models and all of these measurments.Interestingly, the LMC/SMC measurements indicate a kink at L ∼ 3 × 10 4 L ⊙ , with the mass loss rate flattening at smaller luminosities, which does not happen in our models.The observed mass loss rates at low luminosities are closer to the predictions of a Reimer's mass loss prescription.It is likely that a different mass loss mechanism (e.g., mass loss driven by coronal heating due to Alfvén waves, Cranmer and Saar 2011) dominates at low luminosities.
Figure 7 compares our models to the Reimers (Reimers 1975) mass loss rate and we use η = 0.5.Additionally, we show the commonly used empirical mass loss prescription of de Jager et al. (1988), which can be written Finally, we compare to the recent empirical mass loss formula of (Beasor et al. 2020), which can be expressed We do not plot the mass loss prescription of Schröder and Cuntz (2005), which lies between the Reimers and de Jager prescriptions for our models.We can see that the predicted mass loss rates are usually in the range expected for RSGs, 10 −8 M ⊙ /yr ≲ Ṁ ≲ 10 −4 M ⊙ /yr, with canonical mass loss rates of Ṁ ∼ 10 −7 − 10 −5 M ⊙ /yr obtained for models with M ∼ 15 − 20 M ⊙ .We predict strong sensitivity to both luminosity and mass during the RSG phase: low-luminosity ∼ 12 M ⊙ stars are predicted to have very small mass loss rates, Ṁ ∼ 10 −8 − 10 −7 M ⊙ /yr, during core heliumburning, below the Reimer's and Beasor rates.However, when these models expand and brighten during the helium-shell and carbon-burning phases, their mass loss rates increase by more than an order of magnitude to Ṁ > 10 −6 M ⊙ /yr.The strong luminosity dependence of our predicted mass loss rates results from their sensitivity to the ratio of v con /v esc .Brighter RSGS typically have larger luminosities (with slightly larger values of v con ) and larger radii (with slightly smaller values of v esc ), which produces much larger mass loss rates.When RSGs expand and brighten after core helium depletion, their mass loss rates increase substantially, and we predict typical pre-SN mass loss rates of ∼ 10 −5 M ⊙ /yr.For RSGs at the same luminosity/radius, more massive RSGs have much smaller mass loss rates because they have larger values of v esc .This is consistent with the empirical mass loss prescription of Beasor et al. (2020).
Our 15 M ⊙ and 20 M ⊙ models have similar mass loss rates to both the Reimer's and Beasor prescriptions, with typical values of Ṁ ∼ 10 −6 M ⊙ /yr.The models with M ≤ 15 M ⊙ lose a small fraction of their hydrogen envelope before core-collapse, while our 20 M ⊙ model loses about 3 M ⊙ during the RSG phase.The pre-supernova mass loss rate of these models is about ten times larger than during helium burning due to their expansion after helium depletion.
Massive RSGs of M ≳ 30 M ⊙ are predicted to have very large mass loss rates of Ṁ ≳ 10 −5 M ⊙ /yr during core helium-burning, much larger than predicted by the Reimers or Beasor mass loss rates.The 30 M ⊙ model makes it less than halfway through helium burning, living as an RSG for about 1.3 × 10 5 yr before losing most of its hydrogen envelope.The 40 M ⊙ model loses its envelope very quickly, living as an RSG for about 5 × 10 4 yr.In both cases, the models start to evolve into yellow supergiants as they lose mass, where we become less confident in the MLT velocities and our predicted mass loss rates.Those models also require very short time steps, so we terminate their evolution before the end of core heliumburning.It is likely that such stars will lose all of their hydrogen envelopes and turn into Wolf-Rayet stars, but more work will be needed to develop models that reliably follow the stars through that transition.

Implications of mass loss
Our massive models can exhibit exponentially increasing mass loss rate with time.As these stars lose mass, their luminosities/radii remain nearly constant but their masses decrease, decreasing v esc and increasing the mass loss rate, creating a runaway process.Combined with the rarity of such stars due to the IMF, this implies a correspondingly small number of observed very luminous RSGs with large mass loss rates, perhaps consistent with their low observed numbers (Beasor and Smith 2022).In contrast, the Reimers prescription predicts much lower mass loss rates, longer RSG durations, and a larger number of observed very luminous RSGs.The boil-off phenomenon of our models can thus help explain the absence of RSGs observed with luminosities larger than L max ≈ 10 5.5 L ⊙ (Davies et al. 2018), while simultaneously preserving low-mass loss rates of less massive RSGs.Vink and Sabhahit (2023) suggested that the upward kink in mass loss rates at high luminosity (L ≳ 10 4.5 L ⊙ ) is caused by a transition from an optically thin to optically thick wind.This is motivated by an analogy with line-driven winds from hot stars, which show a similar feature as photons begin to scatter multiple times in optically thick winds, allowing them to drive higher mass loss rates.However, as far as we are aware, there is no model of radiatively driven winds in cool stars, e.g., due to photon pressure on molecular absorption lines.While photon pressure on dust helps accelerate RSG winds, the mass loss rate may instead be determined by how much mass can be uplifted to the dust formation radius, so it is not clear that optically thick dust absorption would actually drive higher mass loss rates.Our model can explain the dependence of mass loss on luminosity above the kink, but it cannot explain the dependence below the kink.We speculate that a different mechanism with weaker luminosity dependence dominates mass loss at luminosities below 10 4.5 L ⊙ .Future work will be necessary to distinguish between our suggestion and that of Vink and Sabhahit (2023).
In this work, we do not investigate the effects of metallicity.While our prescription does not explicitly dependent on metallicity, RSG radii (and hence v esc ) do depend on metallicity, which will affect mass loss rates.Preliminary models predict slightly smaller mass loss rates by a factor of less than ∼2 due to smaller radii for LMC metallicity, and very similar to the measurements of (Antoniadis et al. 2024).A competing effect is that lowmetallicity RSGs have lower mass at the same luminosity, increasing mass loss rates.Preliminary SMC metallicity models have similar mass loss rates to the solar metallicity models, given the same luminosity, but lower than the measurements of (Yang et al. 2023).
An apparent prediction of our model is that mass loss does not scale strongly with metallicity.However, there must be some dependence, because very low metallicity stars would not be able to form enough dust to launch the winds at all.A possible prediction of our model is that higher metallicity stars have similar mass loss rates but higher terminal wind speeds due to the larger opacity of dust in their winds.More realistic modeling of dust formation and radiative acceleration will be needed to confirm these hypotheses.

EFFECT ON TYPE II-P SUPERNOVAE
Our model predicts typical densities of ρ ≳ 10 −14 g/cm 3 within three stellar radii of RSGs.This is a few orders of magnitude larger than expected for constant velocity winds with similar mass loss rates (Figure 2).Therefore, we expect signs of CSM interaction in type-II SNe, even for "ordinary" pre-SN mass loss rates of ∼ 10 −6 − 10 −5 M ⊙ /yr.As mentioned above, the total CSM mass of our models is in the range of 10 −3 − 10 −1 M ⊙ , slightly smaller than estimates based on light-curve modeling of type-IIP SNe exhibiting signs of CSM interaction (Morozova et al. 2017;Dessart et al. 2017;Moriya et al. 2018).Our chromospheres qualitatively resemble the less massive but more extended The overlying optical depth profile of the chromosphere after it is fully ionized by supernova shock breakout such that its opacity is κ ≃ 0.35 cm 2 /g.Bottom: Same as top panel, but for time-averaged models.The chromosphere can be optically thick to electron scattering out to more than 2 stellar radii, and its shock breakout radius (where τ = c/v schock ∼ 30) can be greater than 1.2 stellar radii.
CSM structure needed to match both the light curve and narrow emission lines of type II-P (Dessart and Jacobson-Galán 2023), but detailed spectral modeling will be needed to determine whether they can reproduce observations.Here we model the impact of the chromosphere on the light curves of type II-P SNe.
An important property of the CSM is its optical depth.While the star is a RSG, the opacity of this neutral hydrogen/helium is very small (κ ∼ 10 −3 cm 2 /g), and the CSM is optically thin except for large values of v con /v esc .The CSM not greatly obscure the RSG in most cases, unless dust absorption in the wind can redistribute most of the optical/NIR emission to long wavelengths.However, when the star explodes, the shock breakout emission ionizes the CSM, causing its opacity to increase approximately to the Thompson scattering value of κ ≈ 0.35 cm 2 /g.As the CSM recombines (but before being overtaken by the expanding SN ejecta), it may produce narrow flash-ionized emission lines as observed in early spectra of many type II-P SNe.
Figure 8 shows the optical depth of overlying CSM for our models.Assuming Thomson scattering, we predict the CSM to have an optical depth of τ ∼ 10 2 − 10 3 at the star's surface.The shock breakout radius where τ = c/v ej ∼ 30 ranges from ∼ 1.1−1.5 R for our models, which would greatly alter the shock breakout emission and early SN light curve.We do not easily produce CSM with breakout radii larger than 10 14 cm, as observed in ∼1/3 of supernovae according to (Irani et al. 2023).However, -Bolometric and V-band light curves of supernovae from our stellar models with extended chromospheres, given an explosion energy of E = 10 51 erg.For comparison, we also plot the observed light curve of SN 2023ixf (Zimmerman et al. 2023), which showed signs of interaction with circumstellar material, similar to many other type II-P supernovae.
the optical depth of the CSM can be larger than unity out to a few stellar radii (beyond 10 14 cm), which could strongly affect the SN appearance during its first week or two.
We do not perform spectral modeling here, but we note that non-LTE modeling of flash-ionized emission lines has constrained outflow velocities to low values of v wind ≲ 30 km/s in the case of SN 2013cu (e.g., Gräfener and Vink 2016).That event was a hydrogen-poor type-IIb SN for which our models are not directly applicable.In any case, our quasi-hydrostatic chromosphere models predict low outflow/shock speeds within the dust formation radius of exploding type II-P SN progenitors, corresponding to narrow widths of v ≲ 30 km/s of flashionized emission lines, although the widths may be increased by acceleration due to radiation pressure from the subsequent SN (Tsuna et al. 2023).

Light curve modeling
To estimate the effect of the CSM on type-IIP SNe, we perform light curve modeling with the SNEC code (Morozova et al. 2015).The initial models are constructed with our CSM models for v s /v esc = 0.08-0.24,and our time-averaged models with v con /v esc = 0.06-0.10.These are attached to a 15 M ⊙ MESA model at corecollapse, as shown in Figures 2 and 4. We excise the innermost 1.62 M ⊙ of the star corresponding to the oxygen-silicon interface, as suggested by Morozova et al. (2018).The explosions are all simulated with 3000 grid cells, with energy deposited in the innermost 0.1 M ⊙ as a thermal bomb to obtain a final explosion energy of 10 51 erg.The mass of 56 Ni is fixed to 0.05 M ⊙ , but this value does not affect the early part of the light curves, which is our focus.
When the SN shock propagates an extended dense CSM (e.g., the case for v s /v esc = 0.24), the shock downstream is compressed into a thin dense shell due to radiative cooling.By default, SNEC obtains the "observed luminosity" by measuring the local luminosity at the photospheric radius where τ = 2/3.However, due to the thin-shell nature of the shocked region, we find numerical instabilities in the luminosity after the photosphere enters the shocked region (e.g., Fig 4 of Tsuna et al. 2023).To avoid this, we instead measure the bolometric luminosity at the outermost part of the computational region (≈ 3 × 10 16 cm) where photons are free-streaming (τ ≪ 2/3).We use this luminosity and the photospheric radius, and adopt the same bolometric corrections as done in SNEC to obtain the V-band light curves.
Figure 9 shows the bolometric and V-band SN light curves of our models.As the ratio of v s /v esc or v con /v esc increases, the CSM mass increases, and the initial peak in the light curve becomes brighter and longer lasting.This peak arises primarily from shock-cooling emission of the extended CSM.For comparison, we also plot the observed light curve of SN 2023ixf (Zimmerman et al. 2023), which shows prominent signs of interaction with a dense and confined CSM.Models without CSM or with low-density CSM (e.g., v s /v esc = 0.08−0.16)do not come close to matching the observed light curve.Models with denser CSM (e.g., v s /v esc = 0.24 or v con = 0.08 − 0.10) come closer to matching the data, exhibiting broad early peaks in their light curves that can last for more than a week.The agreement is not perfect, but illustrates the point that realistic chromospheres can significantly alter the shape of the early light curve in a manner resembling observations.
Important effects not included in these models are asphericity and intermittency expected for realistic RSG chromospheres.Our light curve models use angleaveraged or time-averaged chrmospheres, whereas we expect the CSM to be very aspherical in real stars, and to be fluctuating greatly in time.This will increase the effect of the CSM in a fraction of events where an energetic convective plume is on the observable side of the star at the time of explosion.In other cases where there is not an energetic convective plume at the time of the explosion, the effect of the CSM will be minor.Additionally, the CSM will be clumpy, which will shorten the diffusion time for some photons and lengthen in for others, broadening the early light curve peak due (Goldberg et al. 2022b).Future multi-dimensional light curve modeling of SNe with clumpy CSM will be needed for detailed comparison with observations.
We note that the model V-band light curves can be affected by SNEC's assumption of gas-radiation equilibrium and a fully thermalized radiation field.We may expect departures from thermal spectra for the luminosity generated around shock breakout, and that powered by prolonged CSM interaction (e.g., Nakar and Sari 2010;Svirski et al. 2012;Morozova et al. 2018;Haynie and Piro 2021;Tsuna et al. 2021).On the other hand, for the shock cooling phase after breakout the thermalization is expected to be more complete (Morozova et al. 2018).Additionally, the morphology of the light curve around the peak is possibly affected by the relatively low resolution within the breakout layer (τ ≈ several tens), most notably for the ones with densest CSM.

Comparison with red supergiant observations
Our model predicts that chromospheric material is clumpy, stochastically varying with time, and moving both inwards and outwards with typical velocities of ∼10 of km/s.These speeds are slightly smaller than observed for gas within 2 stellar radii of Antares and Betel-geuse, moving at velocities of ∼ 20 − 50 km/s (Ohnaka et al. 2017;López Ariste et al. 2022).We also expect that RSGs with larger atmospheric motions should have larger mass loss rates, as tentatively observed (Josselin and Plez 2007).Such observations are typically made though near-infrared emission lines from molecules, which can extend out to a few stellar radii (Arroyo-Torres et al. 2015).It is not clear if our models can account for molecular emission at such large radii, and this should be investigated in future work.
Observations of AGB stars may also be useful, since the same processes are likely to be at work in those stars.For instance, CO emission lines measured with ALMA can be used to infer the density, temperature, and velocity profiles of material in the chromosphere/wind.(Khouri et al. 2024) finds density profiles for R Dor that are very similar to that shown in Figure 2, characterized by a steep drop in density by a factor of ∼ 10 3 within 1.6 R, followed by a more gradual decline in density at larger radii.However, the temperature distribution they infer is somewhat steeper than that of radiative equilibrium (equation 10).More sophisticated thermodynamic modeling and an extension of our model to AGB stars (which are significantly cooler than RSGs) will help refine its predictions.
Our model predicts that most of RSG mass loss occurs through episodic ejection events driven by shocks.These mass ejection events would be accompanied by dust formation events, which may have caused the great dimming of Betelgeuse in 2020 (Montargès et al. 2021;Taniguchi et al. 2022).That event was preceded by radial velocity and molecular absorption disturbances in the photosphere several months before the dust formation (Dupree et al. 2022;Jadlovský et al. 2023).Additional tomography of Betelgeuse showed two strong outgoing shock waves around the time of the great dimming event (Kravchenko et al. 2021;Jadlovský et al. 2024).These observations are qualitatively consistent with our picture of the chromosphere expanding and forming dust during mass ejection events, on a time scale of ∼ R d /v esc ∼ 1 yr.We discuss this further in Section 5.4.

Comparison with prior models
Ours is not the first paper to investigate many of these ideas.Soker (2021) proposed RSGs contain extended "effervescent zones" where some material falls back towards the star while other material is accelerated outwards by photon pressure, arguing this region will be denser than expected from simple wind models.Given the low opacity of chromospheric material, we believe that hydrodynamic propulsion will be more important than photon momentum, at least for material below the dust formation radius.Moriya et al. (2018) examined how the acceleration of RSG winds can affect the CSM density near the stellar surface.They used parameterized β velocity profiles for the wind speed, and then calculated the corresponding wind density profile.If the wind accelerates slowly, the CSM density near the star is correspondingly larger, helping to explain SNe observations requiring dense CSM, in a manner similar to our work.While we find qualitatively similar density profiles to those of a β-law model, the density profiles from our model are likely to be more realistic because they are based on a physical mechanism rather than an arbitrarily chosen velocity profile.These differences may be important when predicting the early evolution of light curves and spectra due to CSM interaction.
It is instructive to compare our model to the recent work of Kee et al. (2021).They consider turbulent pressure support of RSG atmospheres, assuming a constant turbulent speed v turb within the chromosphere.However, they do not consider time time-variable nature of turbulent motions, so their predictions are similar to isothermal density profile of equation 5.In their models, v turb is a free parameter that must be set to ≈ 20 km/s in order to achieve sensible mass loss rates.The problem is that these mass loss rates are incredibly sensitive to the chosen value of v turb : changing v turb from 10 km/s to 20 km/s increases their mass loss rate by about six orders of magnitude (their Figure 3).It seems unlikely that stars naturally fine-tune their turbulent velocities to such an extent.Our model builds on theirs, showing that timeaveraging reduces the sensitivity to v turb , and we link this value to the near-surface convective velocities rather than letting it be a free parameter.
Although sometimes forgotten by the massive star community, models for shock-supported stellar chromospheres have been developed for decades (e.g., Hill and Willson 1979;Bertschinger and Chevalier 1985;Bowen 1988), primarily in the context of radially pulsating AGB stars and Mira variables.Our model is largely a reapplication of previous models to RSGs, but here the shocks are launched by convective motions rather than largescale pulsation.This allows for analytic estimates of chromospheric densities and stellar mass loss rates, which can be readily compared with data or implemented into stellar evolution codes.The most important new aspect of our model is the inclusion of a spectrum of shock speeds v s rather than assuming a single value.The highvelocity tail of this spectrum greatly increases the timeaveraged chromospheric density and mass loss rate.Hill and Willson (1979) found approximately periodic and hydrostatic (in a time-averaged sense) chromospheric solutions as long as the dimensionless ratio P 0 /P (their equation 2) is small.In the relevant limit v con ≪ v esc , this requires As discussed in Section 2, for a a typical period between shocks of P s ∼ 0.5R/v con , we obtain P 0 /P ∼ 8(r/R)(v con /v esc ) 2 , which is less than unity below the dust formation radius for all our RSG models.Hence we expect our assumption of an approximately hydrostatic to be valid except when fairly strong shocks propagate through the chromosphere, which should be investigated more thoroughly with numerical simulations.

Comparison with red supergiant simulations
It is useful to compare our model with simulations of RSGs, such as the 3D radiation hydrodynamical simulations of Goldberg et al. (2022a) (see also Chiavassa et al. 2011a,b;Ma et al. 2023;Chiavassa et al. 2024).These simulations show RMS near-surface convective velocities of ∼10 km/s, similar to those predicted by mixing length theory, for values of α MLT ∼ 2 − 4.These help justify our use MLT to compute the values of v con in our models.However, our predicted mass loss rates are very sensitive to the convective velocities (see Section 5.6), so it will be necessary to directly compute mass loss in simulations in order to test and calibrate our model.
However existing simulations (e.g., those of Goldberg et al. 2022a) may have difficulty correctly predicting the chromospheric density structure.For numerical reasons, the simulations have a density floor, so this material surrounding the star falls onto it.There is also a numerical transient at the start of the simulations that ejects some mass, and some of this mass falls back onto the star over the course of the simulations.These effects cause an accretion shock to form between the stellar surface and the infalling mass (visible in their Figures 6, 10, and 12), which compresses the chromosphere.A similar phenomenon may be present in the CO5BOLD simulations of Wedemeyer et al. (2017) (see their Figure 2).This may compress the chromospheric structure, decreasing the density of material at r ≳ 1.5R.An additional issue is that Goldberg et al. (2022a) assume fully ionized gas in their energy equation and mean molecular weight, whereas real RSGs have primarily neutral atomic H/He in their photospheric and chromosopheric layers.Including recombination energy may be important for the dynamics of material near and above the photosphere.Future simulations that address these issues can better determine the chromospheric structure of RSGs.Some very useful existing simulations are those of AGB stars and their dust-driven winds (Freytag and Höfner 2008;Freytag et al. 2017;Freytag and Höfner 2023).These simulations indeed exhibit RMS velocities of material in and above the photosphere that are nearly constant radius, with RMS radial velocities of ∼5-10 km/s and Mach numbers of order M ∼ 2 − 3. The density profiles of those models also resemble ours in Figure 2, with a sharp drop in density by a few orders of magnitude just above the surface of the star, followed by a more gradual decline at larger radii.Their AGB models have higher values of v con /v esc , cooler surface temperatures, and smaller dust formation radii than we expect for RSGs, so a detailed comparison should await simulations in the RSG regime.

Red supergiant mass loss and variability
Our model predicts that both the chromospheric structure and the mass loss rate of RSGs are dynamic, varying stochastically on a convective turnover time of ∼months.While the plots in this paper are angle-averaged and time-averaged quantities, the actual chromospheric density is expected to exhibit large fluctuations around this mean.Variations of dust spectral features from the wind acceleration region (which are used for mass loss measurements) are likely slower, varying on a time scale of ∼months-years.This could cause inferred mass loss rates to change greatly with time, which could account for differences between different groups' measurements, such as the different mass loss rates of Beasor and Smith (2022) compared to those of van Loon et al. (2005).Future observations should investigate time variations in RSG spectra and inferred mass loss rates.
We can also estimate the likelihood of observing a star with an enhanced CSM density and mass loss rate.As discussed above, the mass loss is dominated by uncom-mon convective episodes with v ∼ √ v con v esc ∼ 3v con .
The probability of a given convective element reaching such velocities is p ∼ e −vesc/2vcon ∼ 10 −3 .For a star with ∼10 surface convective elements and a convective turnover time of ∼100-300 days, we expect to wait ∼30-100 years between major mass ejection events.Each mass loss event would last for a convective turnover time of ∼100 days, ejecting mass at a rate ∼ 10 −4 −10 −3 M ⊙ /yr (see green line in Figure 4), ejecting a total mass of M ej ∼ 3 × 10 −5 − 3 × 10 −4 M ⊙ and about 3×10 −7 −3×10 −6 M ⊙ of dust.For a dust-laden gas opacity of κ d ∼ 3 cm 2 /g, this mass can have an optical depth τ ∼ 1 − 10 at the dust formation radius, temporarily obscuring the star and creating a dimming event similar to that of Betelgeuse in 2020.Such stars would be expected to exhibit great dimming events every century or so, possibly in line with historical evolution of Betelgeuse.Stars with large mass loss rates of 10 −4 M ⊙ /yr (as inferred for the progenitor of SN 2023ixf) would suffer a mass ejection event every few years.Since the gas would remain near the star for a time scale of R d /v wind ∼ 2 yr, the probability of observing such a star with a very dense CSM at the moment of a SN explosion is near unity.

Comparison with supernovae and SN 2023ixf
Observations of type-IIP supernovae have frequently indicated the presence of dense CSM within ∼ 10 15 cm (∼20 stellar radii) of exploding RSGs.This material is revealed by flash-ionized emission lines (e.g., Khazov et al. 2016) from recombination following ionization by shock breakout and continued irradation from the shocked CSM.The emission lines typically fade within ∼5 days after shock-breakout (Bruch et al. 2023) as this CSM is swept up by the SN ejecta that travels outwards at ∼ 10 9 cm/s (∼ 10 14 cm/day).Early light curves of type-IIP SNe often cannot be explained by bare RSGs, but can be explained by interaction with dense CSM.This material has frequently been interpreted and modeled (e.g., Fuller 2017;Morozova et al. 2020;Bruch et al. 2023) as arising from intense pre-supernova mass loss in the final decade of the star's life.
Our model naturally predicts dense chromospheric material within ∼10 stellar radii (∼ 5 × 10 14 cm) of exploding RSGs.This material will be swept up by the SN ejecta in ∼5 days, in line with the disappearance of flashionized emission lines.We believe such material is essentially the chromosphere of the RSG, rather than being mass ejected from the star just before explosion.Pre-SN outbursts likely occur in rarer cases, e.g., SN 2020tlf (Jacobson-Galán et al. 2022), and in type-IIn SNe (Ofek et al. 2014;Strotjohann et al. 2021).Importantly, our model produces dense CSM with only moderate pre-SN mass loss rates of Ṁ ∼ 10 −5 M ⊙ /yr.This is not enough to completely obscure the star (see Section 4), avoiding the problem of progenitor obscuration by superwinds (Davies et al. 2022).
Observations of SN 2023ixf and other type II-P SNe are qualitatively consistent with our model.In SN 2023ixf, modeling of flash-ionized lines (Zhang et al. 2023) suggests CSM densities of ∼ 10 −14 g/cm 3 at ∼5 stellar radii, only slightly above the prediction of our time-averaged model in Figure 4.The outer "edge" of this dense CSM has been inferred to be at ∼10 stellar radii, near the dust formation radius where we expect the CSM density to approach that of a constant-velocity wind model (Bostroem et al. 2023).The average CSM density inferred from that work, ρ ∼ 5 × 10 −14 g/cm 3 , is slightly larger than our time-averaged model, which has a volume-averaged CSM density of ρ ∼ 5 × 10 −15 g/cm 3 between 1.5 stellar radii and the dust formation radius.Zimmerman et al. (2023) infers an effective shock breakout radius for SN 2023ixf of ∼ 2 × 10 14 cm (about 4 stellar radii), substantially larger than our time-averaged model (Figure 8).X-ray observations and SED modeling also suggest large mass loss rates of Ṁ ∼ 10 −4 M ⊙ /yr (Grefenstette et al. 2023;Jencson et al. 2023).We cannot explain this with our models unless the progenitor's mass loss rate at core-collapse rate was larger than average.This is certainly possible since our model predicts stochastic variations in time.Moreover, we expect the CSM to be clumpy such that it is denser and more optically thick on one side of the star, increasing the inferred shock-breakout radius for many viewing angles.
The CSM density at ∼ 10 15 cm is inferred to have a density of ∼ 4×10 −16 g/cm 3 (Zimmerman et al. 2023), which is slightly larger than predicted by our time-averaged model, and again could be reconciled by a model with a larger instantaneous mass loss rate.Zimmerman et al. (2023) interprets this as evidence for a sharp drop in CSM density somewhere above the shock-breakout radius.However, it is not clear from the data that the CSM has a sharp edge: what is required is that the CSM density fall steeper than r −2 from ∼ 10 14 − 10 15 cm, but it could be a smooth drop as predicted by our models.The CSM densities measured from X-ray emission of SN 2023ixf (Grefenstette et al. 2023) are similar to those above, and once again a factor of several larger than our model in Figure 4.
Multiple efforts to model the light curve of SN 2023ixf (Jacobson-Galán et al. 2023;Hiramatsu et al. 2023;Li et al. 2023) have concluded that its early light curve requires dense CSM.For instance, (Jacobson-Galán et al. 2022) infers a density of ρ ∼ 10 −12 cm at ∼2 stellar radii, slightly larger than our time-averaged model in Figure 2.These high densities are often interpreted as evidence for pre-supernova mass loss rates of ∼ 10 −2 − 10 −1 M ⊙ /yr.Similar inferences are often made from flash-ionized emission lines, and have been made for many other type II-P SNe as well (e.g., Yaron et al. 2017;Moriya et al. 2018).While the above CSM density estimates may be reasonable, we believe the extreme mass loss rates are greatly over-estimated.These mass loss rates are obtained by assuming the CSM is moving out at 10s of km/s.Instead, RSG chromospheres within several stellar radii could be approximately hydrostatic and bound, with mean outflow velocities (Figure 5) much slower than assumed by supernova modelers.The actual mass loss rates are correspondingly smaller.Supernova observers and modelers should stick to quoting CSM masses or densities rather than attempting to infer mass loss rates.
Another expectation of our model is that the CSM material is asymmetric and clumpy, stochastically varying in time.This is in line with observations of SN 2023ixf, where the early (and changing) continuum polarization indicates an aspherical CSM (Vasylyev et al. 2023).Spectral modeling of this event also suggests an aspherical CSM (Smith et al. 2023).We suspect an aspherical CSM will also lengthen the duration of shock breakout emission (Goldberg et al. 2022b) and rise to Vband maximum, which may help explain the slow rise of events like SN 2023ixf relative to models (e.g., Figure 9 and Hosseinzadeh et al. 2023).We note that type IIb SN 2013cu was inferred to have a smooth and spherical wind by modeling its narrow emission lines (Gräfener and Vink 2016), but the hydrogen-poor progenitor of that event may have a different mass loss mechanism than our model for hydrogen-rich RSGs.Future hydrodynamic and radiative transfer modelling of our scenario will be needed for detailed predictions.
An interesting and unusual feature of SN 2023ixf is that its progenitor was clearly pulsating at large amplitude (Kilpatrick et al. 2023;Jencson et al. 2023).Many RSGs in the Milky Way and M31 are also observed to pulsate fairly coherently, though often at somewhat smaller amplitudes (Soraisam et al. 2018).Similar pulsations are thought to be crucial to the mass loss of AGB stars by driving shocks through the atmosphere (Höfner and Olofsson 2018;Freytag et al. 2017).Our model could possibly account for this by replacing the convective velocity v con with the pulsational velocity near the photosphere, which can be supersonic and hence larger than v con , driving higher mass loss rates.It is thus possible that the pulsations of SN 2023ixf's progenitor created its unusually large CSM density and mass loss rate, as has been suggested previously (Yoon and Cantiello 2010).

Caveats and uncertainties
The largest uncertainty of our model arises from its sensitive dependence to the convective velocity distribution.Since the chromospheric density is dominated by convective episodes with v ≳ 2v con the density and mass loss from this mechanism depend very sensitively on the high-energy tail of the convective distribution.We have assumed this tail has a Gaussian distribution, as expected for radial velocity variations in isotropic Kolmogorov turbulence.However, the vigorous near-surface convection of RSGs could behave differently.
One difference between these regimes is the convective Mach numbers, which are small in stellar interiors but are of order unity near the surface of a RSG.Convection is neither in the low-Mach number, incompressible regime of Kolmogorov turbulence, nor in the high-Mach number and highly compressible regime of Burgers turbulence.Three-dimensional simulations (e.g., Rabatin and Collins 2023a,b of such turbulence indicate that the convective speeds are fairly well approximated by a Maxwellian distribution, while the density is well approximated by a lognormal distribution, without strong correlations between density and velocity.However, this could be changed by effects of intermittency (Hopkins 2013).Additionally, RSG atmospheres are not isotropic and can reach near-Eddington luminosities, greatly changing the nature of convection (e.g., Goldberg et al. 2022a) and correlations between density and velocity.Further numerical work on non-isotropic, moderate-Mach number, convectively driven turbulence is needed to address this uncertainty.
We neglected radiative forces within the chromosphere, justified by the fact that typical opacities at these temperatures are of order κ ∼ 10 −4 − 10 −2 cm 2 /g (Plez et al. 1992).The corresponding Eddington factor, i.e., the ratio between radiative acceleration and gravity, is f ed = g/g rad = L/L ed = Lκ/(4πGM c) ∼ 10 −4 − 10 −2 at these opacities.This is approximately consistent with radiative transfer modeling of red giant atmospheres (Ireland et al. 2008) that find only small radiative accelerations within the chromosphere, although it has been suggested turbulent motions may decrease line shadowing and increase the effective opacity enough to have a substantial effect (Josselin and Plez 2007).Since our predicted mass loss rates are strongly dependent on v esc and hence on the net effect of gravity and radiation forces, even small radiative forces could substantially increase chromospheric densities and associated mass loss rates, and should be incorporated in more sophisticated versions of our model.
We also neglected the effects of rotation and magnetic fields in our models, which is reasonable given the small rotation rates of RSGs, and the lack of strong surface fields.The cool chromospheric material is almost completely neutral, diminishing the importance of magnetic fields.Furthermore, the thermal and turbulent pressures of our time-averaged model shown in Figure 2 are larger than the magnetic pressure of a dipole field for a surface field strength of 1 G, so fairly large magnetic fields are required in order to dominate the dynamics.Thirumalai andHeyl (2010, 2012) investigated a hybrid magnetic Weber-Davis and dust-driven wind model for red giants and RSGs and its affect on the terminal wind speeds.However, the mass loss rate is imposed in such models (rather than being calculated from first principles), so such models do not solve the wind launching problem below the dust formation radius.
Finally, our model for the dust-driven wind above the dust formation radius is undoubtedly too simple.More realistic dust chemistry, grain formation, and radiative driving should be developed for a better model of the wind material.Additionally, resonant drag instabilities will occur in the dust-driven ouflow (Squire and Hopkins 2018;Moseley et al. 2019;Steinwandel et al. 2022), which may affect the outflow speed and appearance.Understanding such phenomena may be needed not only for mass loss models, but also for more accurate observational mass loss rate estimates.

CONCLUSIONS
We have developed a simple model for material supported by shock waves in the chromospheres of red supergiants (RSGs).These chromospheres extend roughly from the stellar surface to the dust formation radius at ∼5 stellar radii.Gravity is balanced by momentum deposition from outgoing shock waves, whose effect is similar to turbulent pressure support.For our assumption of nearly constant shock speeds (justified by simulation results), the density profile is similar to that of an isothermal atmosphere with scale height v 2 con /g ≪ R, which would cause a steeply decreasing density for shock speeds comparable to the surface convective speeds v con .However, we also showed that averaging over a stochastically varying convective velocity spectrum greatly increases the chromospheric density, causing it to fall off less steeply with radius.
Using the time-averaged chromospheric density at the dust formation radius, we construct a simple and analytic estimate for the mass loss rate of RSGs (equation 15).This yields mass loss rates approximately compatible with observations of luminous RSGs (Figure 7).The predicted mass loss rates scale exponentially with the ratio of v con /v esc , resulting in larger mass loss rates for larger radii, smaller mass, or higher convective velocities.These models can potentially explain: i) mass loss rates of ∼ 10 −6 M ⊙ /yr for 15 − 20 M ⊙ RSGs, ii) substantially lower mass loss rates for less massive RSGs as found by recent observations, iii) the absence of RSGs for luminosities L ≳ 3 × 10 5 L ⊙ due to rapid boil-off of their hydrogen envelopes.Our models also predict presupernova mass loss rates of ∼ 4×10 −7 −4×10 −5 M ⊙ /yr, a factor 10-100 times larger than mass loss rates during core helium burning.
Our models predict chromospheric masses of ∼ 10 −2 M ⊙ , most of which is confined within one stellar radius above the photosphere.Nonetheless, the typical chromospheric densities of our models can be a few orders of magnitude larger than predicted by a constantvelocity wind model.This chromospheric material likely helps explain i) flash-ionized emission lines observed in the first week of many type II SNe ii) early peaks in the light curves of these SNe that cannot be explained for bare RSG models iii) extended shock breakout emission (and larger shock breakout radius) of type II SNe.
Because convective turbulence varies stochastically, we expect RSGs chromospheres to vary stochastically as well, meaning they are likely to be spatially asymmetric at any moment in time, and their total mass varies greatly in time.This may account for the time-varying molecular emission that is interferometrically observed above the surfaces of RSGs such as Betelgeuse and Antares.It may also explain the great dimming event of Betelgeuse as a moment of enhanced convective speeds, leading to enhanced mass loss, dust formation, and obscuration.Such stochasticity is also expected for type II SNe, meaning that the CSM densities (and flash ionized line strengths, etc.) will likely vary greatly from one SN to the next, even for progenitors of the same mass.
Our chromospheric model is extremely simple and will need to be refined and improved for more robust predictions.Since the chromospheric mass scales exponentially with the convective velocity v con , the most important uncertainty of our model is the actual velocity spectrum of convective motions and their radial dependence within the chromosphere.The model will also be improved with a better treatment of dust formation and radiative driving of the wind around the dust formation radius.These uncertainties can be addressed with future simulations of RSG chromospheres, continued observations of RSGs and their mass loss rates, and early observations of type II SNe.This paper was built using the Open Journal of Astrophysics L A T E X template.The OJA is a journal which provides fast and easy peer review for new papers in the astro-ph section of the arXiv, making the reviewing process simpler for authors and referees alike.Learn more at http://astro.theoj.org.

Ṁ
Fig. 2.-Top: Density profiles of a chromosphere model (solid lines) for representative values of vs/vesc for a 15 M ⊙ MESA model of a red supergiant just before core-collapse.This model has an RMS convective velocity of vcon/vesc ≃ 0.08, and a mass loss rate of Ṁ ≈ 10 −5 M ⊙ /yr.Dashed lines show density profiles for winds of the indicated mass loss rates and a constant wind speed of v wind = 30 km/s.Bottom: Time-averaged density profiles of the same model, for different values of the RMS convective velocity that are representative of red supergiants.

Fig. 4 .
Fig. 4.-Top: Same as Figure 2, but now showing the transition from chromosphere to wind at a dust formation radius of R d = 5R (dashed black line).Near the stellar surface, the density falls roughly exponentially, and it falls roughly as r −2 at larger radii.Bottom: Same as top panel, but showing the time-averaged density profiles.
Fig.7.-Mass loss rates for the same models as Figure7.The black symbols are the Reimer's mass loss rates of the same models for η = 0.5, while the dark gray symbols are for thede Jager et al. 1988 prescription, and dark gray symbols are the empirical prescription ofBeasor et al. 2020.Our model predicts a steeper scaling of mass loss with luminosity than any of these prescriptions.
Fig.8.-Top:The overlying optical depth profile of the chromosphere after it is fully ionized by supernova shock breakout such that its opacity is κ ≃ 0.35 cm 2 /g.Bottom: Same as top panel, but for time-averaged models.The chromosphere can be optically thick to electron scattering out to more than 2 stellar radii, and its shock breakout radius (where τ = c/v schock ∼ 30) can be greater than 1.2 stellar radii.
Fig.9.-Bolometric and V-band light curves of supernovae from our stellar models with extended chromospheres, given an explosion energy of E = 10 51 erg.For comparison, we also plot the observed light curve of SN 2023ixf(Zimmerman et al. 2023), which showed signs of interaction with circumstellar material, similar to many other type II-P supernovae.