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

Jim Fuller1 and Daichi Tsuna1,2
1TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA
91125, USA
2Research Center for the Early Universe (RESCEU), School of Science, The
University of Tokyo, Bunkyo, Tokyo 113-0033, Japan Email<EMAIL_ADDRESS>
###### Abstract

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 may 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.

###### keywords:

##  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_{\odot}/{\rm yr}\lesssim\dot{M}\lesssim
10^{-4}M_{\odot}/{\rm 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_{\odot}\lesssim M\lesssim 40\,M_{\odot}$) during their post-main
sequence evolution, with large radii of $300\,R_{\odot}\lesssim R\lesssim
1000\,R_{\odot}$, high luminosities of $10^{4}\,L_{\odot}\lesssim L\lesssim
3\times 10^{5}\,L_{\odot}$, and cool surface temperatures of $3500\,{\rm
K}\lesssim T_{\rm eff}\lesssim 4500\,{\rm K}$ (e.g., Levesque et al. 2005).
Their low escape speeds $v_{\rm esc}\sim 100\,{\rm km}/{\rm 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 & Wickramasinghe 1962), which forms well above their
photospheres where the temperature drops to $T\lesssim 1500\,{\rm K}$ (Field
1974; Höfner & 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_{\rm d}\gtrsim
5\,R$, and the mechanisms by which material is lofted up to the dust formation
radius are highly uncertain. This wind-launching 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 & 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 & Saar 2011) so that its
temperature maintains radiative equilibrium with the star and decreases
approximately as $T\propto r^{-1/2}$ in the chromosphere. Evermore
sophisticated radiation hydrodynamics simulations of this process (Freytag &
Höfner, 2008; Freytag et al., 2017; Freytag & 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 surroundings (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 $\sim\\!10^{9}\,{\rm
cm}/{\rm s}$, the outgoing SN ejecta takes $\sim$5 days to expand from the
star’s photosphere out to the dust formation radius of $R_{\rm d}\\!\sim
10\,R\sim\\!5\times 10^{14}\,{\rm 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., 2021, 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,
$\dot{M}\sim 10^{-6}\,M_{\odot}/{\rm yr}$, $v_{\rm wind}\sim 30\,{\rm km}/{\rm
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 weeks-
months before explosion (Woosley & 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.

##  2 Chromospheric Dynamics

Figure 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 $\sim$5 stellar radii forms dust, increasing its opacity and creating a
radiation-driven wind.

###  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_{\rm con}\sim 10\,{\rm km}/{\rm s}$ (Goldberg et al.,
2022a). In comparison, the escape speeds from the are

$\displaystyle v_{\rm esc}$ $\displaystyle=\sqrt{2GM/R}$ $\displaystyle\simeq
100\,{\rm
km/s}\bigg{(}\frac{M}{14M_{\odot}}\bigg{)}^{\\!1/2}\bigg{(}\frac{R}{500R_{\odot}}\bigg{)}^{\\!-1/2}\,,$
(1)

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 the photosphere, $\Delta R\sim R(v_{\rm con}/v_{\rm
esc})^{2}\sim R/100$, so the density must fall sharply.

Simulations show that turbulent pressure support dominates over thermal
pressure support for material above the photosphere. For a turbulent pressure
$P_{\rm turb}\sim\rho v_{\rm turb}^{2}$, where $\rho$ is the density and
$v_{\rm turb}$ is the RMS turbulent velocity, we naively expect a
chromospheric scale height of $H\sim P_{\rm turb}/(\rho g)\sim v_{\rm
con}^{2}/g\sim 2r(v_{\rm turb}/v_{\rm esc})^{2}$. This is a factor of $\sim$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

$\frac{\partial}{\partial
t}\big{(}\rho{\vec{v}}\big{)}+\nabla\cdot\big{(}\rho{\vec{v}}{\vec{v}}\big{)}=-\nabla
P+\rho{\vec{g}}\,.$ (2)

The $\rho{\vec{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 & 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 $\gamma\simeq 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_{\rm s}$ (e.g., Shu 1992). Hence, the radial momentum
flux of an outward going shock wave is $\rho v_{\rm s}^{2}$, where $v_{s}$ is
the shock speed.

Applying these simplifications to equation 2, hydrostatic equilbrium requires

$\frac{1}{r^{2}}\frac{d}{dr}\big{(}\rho r^{2}v_{\rm
s}^{2}\big{)}=-\frac{GM\rho}{r^{2}}\,$ (3)

where $r$ is the radial coordinate. Here, $\rho$ 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 $\mathcal{M}\sim 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 & McKee
(1999), a shock wave propagating through a star accelerates the material to a
velocity of

$v_{s}\sim
v_{0}\bigg{(}\frac{m}{m_{0}}\bigg{)}^{\\!\beta-1/2}\bigg{(}\frac{\rho}{\rho_{0}}\bigg{)}^{\\!-\beta}\bigg{(}\frac{r}{r_{0}}\bigg{)}^{\\!-3\beta}\,.$
(4)

Here, $v_{0}$, $m_{0}$, $\rho_{0}$, and $r_{0}$ are the shock velocity, mass
coordinate, density, and radius at an arbitrary reference location. The value
of $\beta$ depends on the equation of state of the gas (Sakurai, 1960), but is
approximately $\beta\approx[2+\sqrt{2\gamma/(\gamma-1)}]^{-1}$ for strong
shocks (Whitham, 2011), which gives $\mu\simeq 0.2$ for $\gamma=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, $\gamma\simeq 1$ and $\beta\simeq 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}\sim 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 & 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_{\rm con}\sim 0.5R/v_{\rm
con}$ such that the separation between shocks is $\sim 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
$\dot{p}=d(4\pi r^{2}\rho v_{\rm s})/dt$. For a time-steady state, the time
derivative is $d/dt\rightarrow v_{s}d/dr$. If $v_{\rm s}$ is constant as a
function of radius, it can be moved inside the derivative. The momentum
deposition per unit volume is thus $\dot{p}/(4\pi r^{2})=(1/r^{2})d(\rho
r^{2}v_{\rm s}^{2})/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_{\rm s}$.

###  2.2 Chromospheric Structure

Proceeding with the assumption that $v_{\rm s}$ is constant as a function of
radius, equation 3 can be solved to yield the density profile

$\rho=\rho(R)\bigg{(}\frac{R}{r}\bigg{)}^{2}e^{-\frac{v_{\rm esc}^{2}}{2v_{\rm
s}^{2}}\big{(}1-\frac{R}{r}\big{)}}\,.$ (5)

This is similar to the structure of an isothermal atmosphere, but with the
sound speed replaced by $v_{\rm 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 & 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_{\rm 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

$p(v_{\rm s})=\sqrt{\frac{2}{\pi}}\frac{1}{v_{\rm con}}e^{-\frac{v_{\rm
s}^{2}}{2v_{\rm con}^{2}}}\,,$ (6)

where $v_{\rm con}$ is the typical convective velocity. From this definition,
the mean turbulent velocity is $v_{\rm con}$. We use a Gaussian rather than
Maxwellian distrubtion because we consider only the radial component of the
convective velocities.

Figure 2:  Top: Density profiles of a chromosphere model (solid lines) for
representative values of $v_{\rm s}/v_{\rm esc}$ for a $15\,M_{\odot}$ MESA
model of a red supergiant just before core-collapse. This model has an RMS
convective velocity of $v_{\rm con}/v_{\rm esc}\simeq 0.08$, and a mass loss
rate of $\dot{M}\approx 10^{-5}\,M_{\odot}/{\rm yr}$. Dashed lines show
density profiles for winds of the indicated mass loss rates and a constant
wind speed of $v_{\rm wind}=30\,{\rm km}/{\rm s}$. Bottom: Time-averaged
density profiles of the same model, for different values of the RMS convective
velocity that are representative of red supergiants.

Figure 3:  Top: Profile of overlying chromospheric mass as a function of
radius corresponding to the models in Figure 2. We expect chromospheres of
$\sim\\!10^{-2}-10^{-1}\,M_{\odot}$, mostly within $\sim$2 stellar radii.
Bottom: Time-averaged mass profiles of the same model.

If we assume that $v_{\rm s}$ is constant throughout the chromosphere, the
time-averaged density profile is

$\displaystyle\rho_{\rm av}(r)$ $\displaystyle=\int^{\infty}_{0}p(v_{\rm
s})\,\rho(r)\,dv_{\rm s}$
$\displaystyle=\sqrt{\frac{2}{\pi}}\rho(R)\bigg{(}\frac{R}{r}\bigg{)}^{2}\int^{\infty}_{0}e^{-\frac{v_{\rm
s}^{2}}{2v_{\rm con}^{2}}-\frac{v_{\rm esc}^{2}}{2v_{\rm
s}^{2}}\big{(}1-\frac{R}{r}\big{)}}\frac{dv_{\rm s}}{v_{\rm con}}$
$\displaystyle=\rho(R)\bigg{(}\frac{R}{r}\bigg{)}^{2}e^{-\frac{v_{\rm
esc}\sqrt{1-R/r}}{v_{\rm con}}}\,.$ (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_{\rm esc}/v_{\rm
con}\sim 10$ is a large number, so the time-averaged density profile of
equation 2.2 falls off much more slowly than that of equation 5. The integrand
of equation 2.2 is maximized where $v_{\rm s}\sim\sqrt{v_{\rm esc}v_{\rm
con}}\sim 3v_{\rm 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 $\sim$2 sigma larger than their usual values.

Figure 2 shows the density profiles of equation 5 and equation 2.2 above a RSG
photosphere. The underlying model is a $13.8\,M_{\odot}$ RSG with
$R=844R_{\odot}$ generated using MESA (Paxton et al., 2011; Paxton et al.,
2013, 2015). This model had a ZAMS mass of $15\,M_{\odot}$ and has been
evolved to core-collapse. We use MLT with $\alpha_{\rm MLT}=2.5$ to measure
the convective velocity, which we evaluate where the optical depth is
$\tau=P_{\rm rad}c/(Pv_{\rm con})$, where $P_{\rm 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 sound speed $c_{s}$. This model has $v_{\rm con}\simeq
0.08v_{\rm esc}$ near the star’s surface.

We plot the density profiles of equation 5 for $v_{\rm s}/v_{\rm esc}=$0.08,
0.16, and 0.24 to demonstrate how moments of large $v_{\rm s}$ greatly
increase the chromospheric density. We also plot the time-averaged density
profile of equation 2.2 for values of $v_{\rm con}/v_{\rm 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 high-mass RSGS or pre-SN RSGs at the
upper end. For this $15\,M_{\odot}$ model, the corresponding range in mass
loss rate is $\dot{M}\sim 10^{-7}-10^{-4}\,M_{\odot}/{\rm yr}$ (see Section
3), demonstrating its senstivity to $v_{\rm con}/v_{\rm 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 $\dot{M}_{\rm wind}$ has a density

$\rho_{\rm wind}=\frac{\dot{M}_{\rm wind}}{4\pi r^{2}v_{\rm wind}}\,.$ (8)

We plot $\rho_{\rm wind}$ for $v_{\rm 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 $\dot{M}_{\rm wind}=10^{-4}\,M_{\odot}$/yr. Between
1-2 stellar radii, our chromospheres may have densities comparable to an
extreme wind of $\dot{M}_{\rm wind}=10^{-2}\,M_{\odot}$/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

$M_{\rm CSM}(r)=\int^{R_{\rm d}}_{r}4\pi\rho r^{2}dr\,$ (9)

where $R_{\rm 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_{\rm d}$ (see next section) of
$\sim\\!10^{-2}\,M_{\odot}$, with large variation for different values of
$v_{\rm s}/v_{\rm esc}$ or $v_{\rm con}/v_{\rm esc}$. Most of this mass is
near the photosphere, but the time-averaged mass above two stellar radii can
be $\sim\\!10^{-3}\,M_{\odot}$. For comparison, a constant wind model with a
similar mass loss rate of $\dot{M}=10^{-5}\,M_{\odot}$/yr contains only
$\sim\\!10^{-5}\,M_{\odot}$ within two stellar radii.

##  3 Mass Loss

###  3.1 Dust Formation and Wind Acceleration

Figure 4:  Top: Same as Figure 2, but now showing the transition from
chromosphere to wind at a dust formation radius of $R_{\rm 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.  Figure 5:  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.

If chromospheric material extends to large radii, it can cool enough to form
large amounts of dust. Radiation 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)

$\displaystyle T(r)$
$\displaystyle=T_{*}\bigg{[}\frac{1-\big{(}1-(R/r)^{2}\big{)}^{1/2}}{2}\bigg{]}^{\\!1/4}$
$\displaystyle\approx\frac{T_{*}}{2^{1/2}}\bigg{(}\frac{R}{r}\bigg{)}^{\\!1/2}\qquad{\rm
for}\,r\gtrsim 2R\,.$ (10)

where $T_{*}$ is the photospheric temperature. For a given dust formation
temperature $T_{\rm d}$, the corresponding dust formation radius is

$R_{\rm d}\approx\frac{R}{2}\bigg{(}\frac{T_{*}}{T_{\rm d}}\bigg{)}^{\\!2}\,.$
(11)

For RSG surface temperatures of $T_{*}\approx 3600\,{\rm K}$ and dust
formation temperatures of $T_{\rm d}\sim 1200\,{\rm K}$, dust condenses at a
dust formation radius of $R_{\rm d}\sim 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_{\rm d}$ and
falling off as $1/r^{2}$ above the dust formation radius. We also neglect
turbulent pressure support above $R_{\rm d}$ so that the momentum equation in
this region is

$v\frac{dv}{dr}=-\frac{GM}{r^{2}}+\frac{f_{\rm Ed}GM}{r^{2}}\qquad{\rm
for}\,\,r>R_{\rm d}\,.$ (12)

Here, the constant $f_{\rm Ed}=L/L_{\rm Ed}$ determines the magnitude of the
radiation force on the dust-gas mixture and will be computed below. For an
initial velocity $v(R_{\rm d})=v_{\rm d}$ at the dust formation radius, the
solution to this equation is

$v_{\rm out}=\bigg{[}v_{\rm d}^{2}+x_{\rm d}GM\bigg{(}\frac{1}{R_{\rm
d}}-\frac{1}{r}\bigg{)}\bigg{]}^{1/2}\qquad{\rm for}\,\,r>R_{\rm d}$ (13)

where $x_{\rm d}=2(f_{\rm Ed}-1)$. The velocity at $r\rightarrow\infty$ is
then

$v_{\infty}=\sqrt{v_{\rm d}^{2}+x_{\rm d}GM/R_{\rm d}}\,.$ (14)

To match our solution above and below the dust formation radius, we choose an
outflow velocity $v_{\rm d}=v_{\rm 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_{\rm con}$, and the advection term
becomes important in the momentum equation. The mass loss rate is then

$\displaystyle\dot{M}$ $\displaystyle\approx 4\pi R_{d}^{2}\rho(R_{d})\,v_{\rm
con}$ $\displaystyle\approx 4\pi R^{2}\rho(R)v_{\rm con}\,e^{-\frac{v_{\rm
esc}\sqrt{1-R/R_{\rm d}}}{v_{\rm con}}}\,.$ (15)

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_{\rm s}=GM/v_{\rm
con}^{2}\sim 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_{\infty}$ to constrain $x_{\rm
d}$, but one could use dust opacities and mass fractions to calculate $x_{\rm
d}$ and then predict the corresponding value of $v_{\infty}$. From equation
14, we have $x_{\rm d}=(R_{\rm d}/GM)(v_{\infty}^{2}-v_{\rm d}^{2})$. In this
work, we use $v_{\infty}=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 & Josselin 2011). This requires $x_{\rm d}\sim 2$, approximately what
we expect for realistic dust-laden gas opacities of $\kappa_{\rm d}\sim
5\,{\rm cm^{2}}/{\rm g}$ (Höfner et al., 2003).

To compute the time-averaged outflow speed below the dust formation radius, we
simply require $\dot{M}$ to be constant with radius such that $v_{\rm
out}=\dot{M}/(4\pi\rho 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

$\rho(r)=\frac{\dot{M}}{4\pi r^{2}}\bigg{[}v_{\rm d}^{2}+x_{\rm
d}GM\bigg{(}\frac{1}{R_{\rm d}}-\frac{1}{r}\bigg{)}\bigg{]}^{-1/2}\qquad{\rm
for}\,\,r>R_{\rm d}\,.$ (16)

As mentioned above, dust formation begins at temperatures ranging from
$T\simeq 1000-2000\,{\rm K}$, with different species condensing at different
temperatures. Here we set the dust formation radius to $R_{\rm d}=5R$, a
reasonable value for RSGs. Different choices for $R_{\rm d}$ would result in
slightly different density/velocity profiles above $R_{\rm d}$. The mass loss
rate of our model is only weakly sensitive to the value of $R_{\rm d}$ because
the factor of $R_{\rm d}^{-2}$ in the density profile (equation 5) cancels out
the factor of $R_{\rm d}^{2}$ in the mass loss rate (equation 3.1). 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_{\rm con}/v_{\rm esc}$. The
density profiles are the same as Figure 2 below the dust formation radius, but
we now extend the plot above $R_{\rm d}$. The velocity profile has low outflow
velocities below $R_{\rm 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_{\rm d}$, and then it approaches the constant-velocity solution
above $R_{\rm 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_{\rm d}$. This is caused by the upward kink in
the outflow velocity at $R_{\rm 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.

###  3.2 Red Supergiant Mass Loss Rates

Figure 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_{\rm
eff}<4200\,{\rm 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_{\odot}$ and $40\,M_{\odot}$ 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\\!\lesssim 2\times 10^{4}\,L_{\odot}$.  Figure 7:  Mass loss rates for the
same models as Figure 7. The black symbols are the Reimer’s mass loss rates of
the same models for $\eta=0.5$, while the dark gray symbols are for the de
Jager et al. 1988 prescription, and dark gray symbols are the empirical
prescription of Beasor et al. 2020. Our model predicts a steeper scaling of
mass loss with luminosity than any of these prescriptions.

We implement the mass loss rate of the previous section (equation 3.1) into
stellar evolution models generated with MESA as described above. At surface
temperatures above $4500\,{\rm K}$ we use the “Dutch” wind prescription with
scaling factor $\eta=0.5$. At cooler surface temperatures, we use our wind
prescription.

Figure 6 shows the predicted mass loss rates for RSGs with masses of
$10-40\,M_{\odot}$. We only plot the RSG phase where $T_{\rm eff}<4200\,{\rm
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_{\rm 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\times 10^{4}\,L_{\odot}$ is quite similar
between our models and all of these measurments. Interestingly, the LMC/SMC
measurements indicate a kink at $L\sim 3\times 10^{4}\,L_{\odot}$, 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 & Saar 2011) dominates at low luminosities.

Figure 7 compares our models to the Reimers (Reimers, 1975) mass loss rate

$\dot{M}_{\rm Reim}=\eta\,4\times
10^{-13}\frac{(L/L_{\odot})(R/R_{\odot})}{M/M_{\odot}}M_{\odot}/{\rm yr}\,,$
(17)

and we use $\eta=0.5$. Additionally, we show the commonly used empirical mass
loss prescription of de Jager et al. (1988), which can be written

$\dot{M}_{\rm deJa}=4.5\times
10^{-6}\bigg{(}\frac{L}{10^{5}L_{\odot}}\bigg{)}^{\\!1.769}\bigg{(}\frac{T_{\rm
eff}}{4000\,{\rm K}}\bigg{)}^{\\!-1.676}\,M_{\odot}/{\rm yr}$ (18)

Finally, we compare to the recent empirical mass loss formula of (Beasor et
al., 2020), which can be expressed

$\dot{M}_{\rm
Beas}=10^{-3.5}10^{-3M/(20M_{\odot})}\bigg{(}\frac{L}{L_{\odot}}\bigg{)}^{\\!3.6}\,M_{\odot}/{\rm
yr}\,.$ (19)

We do not plot the mass loss prescription of Schröder & 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_{\odot}/{\rm yr}\lesssim\dot{M}\lesssim
10^{-4}\,M_{\odot}/{\rm yr}$, with canonical mass loss rates of $\dot{M}\sim
10^{-7}-10^{-5}\,M_{\odot}/{\rm yr}$ obtained for models with
$M\\!\sim\\!15-20\,M_{\odot}$. We predict strong sensitivity to both
luminosity and mass during the RSG phase: low-luminosity
$\sim\\!12\,M_{\odot}$ stars are predicted to have very small mass loss rates,
$\dot{M}\sim 10^{-8}-10^{-7}\,M_{\odot}/{\rm yr}$, during core helium-burning,
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
$\dot{M}>10^{-6}\,M_{\odot}/{\rm yr}$.

The strong luminosity dependence of our predicted mass loss rates results from
their sensitivity to the ratio of $v_{\rm con}/v_{\rm esc}$. Brighter RSGS
typically have larger luminosities (with slightly larger values of $v_{\rm
con}$) and larger radii (with slightly smaller values of $v_{\rm 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 $\sim 10^{-5}\,M_{\odot}/{\rm yr}$. For RSGs
at the same luminosity/radius, more massive RSGs have much smaller mass loss
rates because they have larger values of $v_{\rm esc}$. This is consistent
with the empirical mass loss prescription of Beasor et al. (2020).

Our $15\,M_{\odot}$ and $20\,M_{\odot}$ models have similar mass loss rates to
both the Reimer’s and Beasor prescriptions, with typical values of
$\dot{M}\sim 10^{-6}\,M_{\odot}/{\rm yr}$. The models with $M\leq
15\,M_{\odot}$ lose a small fraction of their hydrogen envelope before core-
collapse, while our $20\,M_{\odot}$ model loses about $3\,M_{\odot}$ 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\\!\gtrsim\\!30\,M_{\odot}$ are predicted to have very large
mass loss rates of $\dot{M}\gtrsim 10^{-5}\,M_{\odot}/{\rm yr}$ during core
helium-burning, much larger than predicted by the Reimers or Beasor mass loss
rates. The $30\,M_{\odot}$ model makes it less than halfway through helium
burning, living as an RSG for about $1.3\times 10^{5}\,{\rm yr}$ before losing
most of its hydrogen envelope. The $40\,M_{\odot}$ model loses its envelope
very quickly, living as an RSG for about $5\times 10^{4}\,{\rm 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 helium-burning. 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.

###  3.3 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_{\rm 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 & 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_{\rm max}\approx 10^{5.5}\,L_{\odot}$ (Davies
et al., 2018), while simultaneously preserving low-mass loss rates of less
massive RSGs.

Vink & Sabhahit (2023) suggested that the upward kink in mass loss rates at
high luminosity ($L\gtrsim 10^{4.5}\,L_{\odot}$) 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_{\odot}$. Future work will be necessary to
distinguish between our suggestion and that of Vink & 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_{\rm 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 $\sim$2 due to smaller radii for LMC metallicity, and very
similar to the measurements of (Antoniadis et al., 2024). A competing effect
is that low-metallicity 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.

Figure 8:  Top: The overlying optical depth profile of the chromosphere after
it is fully ionized by supernova shock breakout such that its opacity is
$\kappa\simeq 0.35\,{\rm cm}^{2}/{\rm 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 $\tau=c/v_{\rm schock}\sim 30$) can be greater than 1.2 stellar radii.

##  4 Effect on Type II-P Supernovae

Our model predicts typical densities of $\rho\gtrsim 10^{-14}\,{\rm g}/{\rm
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
$\sim\\!10^{-6}-10^{-5}\,M_{\odot}/{\rm yr}$. As mentioned above, the total
CSM mass of our models is in the range of $10^{-3}-10^{-1}\,M_{\odot}$,
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 CSM structure needed to match both the light curve
and narrow emission lines of type II-P (Dessart & 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.

Figure 9:  Bolometric and V-band light curves of supernovae from our stellar
models with extended chromospheres, given an explosion energy of
$E=10^{51}\,{\rm 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.

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 ($\kappa\sim
10^{-3}\,{\rm cm}^{2}/{\rm g})$, and the CSM is optically thin except for
large values of $v_{\rm con}/v_{\rm 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 $\kappa\approx 0.35\,{\rm
cm^{2}}/{\rm 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 $\tau\sim
10^{2}-10^{3}$ at the star’s surface. The shock breakout radius where
$\tau=c/v_{\rm ej}\sim 30$ ranges from $\sim\\!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}\,{\rm cm}$, as observed in $\sim$1/3 of supernovae according to
(Irani et al., 2023). However, the optical depth of the CSM can be larger than
unity out to a few stellar radii (beyond $10^{14}\,{\rm 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_{\rm wind}\lesssim 30\,{\rm km/s}$ in the case of SN 2013cu (e.g.,
Gräfener & 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\lesssim 30\,{\rm km/s}$ of flash-ionized emission lines, although the
widths may be increased by acceleration due to radiation pressure from the
subsequent SN (Tsuna et al., 2023a).

###  4.1 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_{\rm s}/v_{\rm esc}=0.08$—$0.24$, and
our time-averaged models with $v_{\rm con}/v_{\rm esc}=0.06$—$0.10$. These are
attached to a 15 $M_{\odot}$ MESA model at core-collapse, as shown in Figures
2 and 4. We excise the innermost $1.62\ M_{\odot}$ 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_{\odot}$ as a thermal bomb to obtain a final explosion
energy of $10^{51}$ erg. The mass of 56Ni is fixed to $0.05\ M_{\odot}$, 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_{\rm
s}/v_{\rm 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
$\tau=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., 2023b). To avoid this, we
instead measure the bolometric luminosity at the outermost part of the
computational region ($\approx 3\times 10^{16}$ cm) where photons are free-
streaming ($\tau\ll 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_{\rm s}/v_{\rm esc}$ or $v_{\rm con}/v_{\rm 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_{\rm s}/v_{\rm esc}=0.08-0.16$) do not come close to matching the
observed light curve. Models with denser CSM (e.g., $v_{\rm s}/v_{\rm
esc}=0.24$ or $v_{\rm 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 angle-averaged 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 & Sari, 2010; Svirski et al., 2012; Morozova et al., 2018; Haynie
& 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 ($\tau\approx$ several tens), most notably for the ones
with densest CSM.

##  5 Discussion

###  5.1 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 $\sim$10 of km/s. These speeds are slightly smaller than
observed for gas within 2 stellar radii of Antares and Betelgeuse, moving at
velocities of $\sim\\!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 & 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 $\sim\\!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 3.1). 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. 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), qualitatively consistent with our picture of the chromosphere
responding to energetic convection below on a dynamical time scale of
$\sim\\!R_{\rm d}/v_{\rm esc}\sim 1\,{\rm yr}$. We discuss this further in
Section 5.4.

###  5.2 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 $\beta$ 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 $\beta$-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_{\rm 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_{\rm turb}$ is a free parameter that must be set to $\approx
20\,{\rm km}/{\rm 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_{\rm turb}$: changing $v_{\rm 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 time-
averaging reduces the sensitivity to $v_{\rm 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 &
Willson 1979; Bertschinger & 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 large-scale 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
high-velocity tail of this spectrum greatly increases the time-averaged
chromospheric density and mass loss rate.

Hill & 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_{\rm con}\ll
v_{\rm esc}$, this requires

$\frac{P_{0}}{P}\simeq\frac{4rv_{\rm con}}{v_{\rm esc}^{2}P_{\rm s}}\,.$ (20)

As discussed in Section 2, for a a typical period between shocks of $P_{\rm
s}\sim 0.5R/v_{\rm con}$, we obtain $P_{0}/P\sim 8(r/R)(v_{\rm con}/v_{\rm
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.

###  5.3 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. 2011b; Chiavassa et al. 2011a; Ma et al. 2023; Chiavassa et
al. 2024). These simulations show RMS near-surface convective velocities of
$\sim$10 km/s, similar to those predicted by mixing length theory, for values
of $\alpha_{\rm MLT}\sim 2-4$. These help justify our use MLT to compute the
values of $v_{\rm 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\\!\gtrsim 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 & Höfner, 2008; Freytag et al., 2017; Freytag & 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 $\sim$5-10 km/s and Mach numbers of order $\mathcal{M}\sim 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_{\rm con}/v_{\rm 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.

###  5.4 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 $\sim$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 $\sim$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 & 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
uncommon convective episodes with $v\sim\sqrt{v_{\rm con}v_{\rm esc}}\sim
3v_{\rm con}$. The probability of a given convective element reaching such
velocities is $p\sim e^{-v_{\rm esc}/2v_{\rm con}}\sim 10^{-3}$. For a star
with $\sim$10 surface convective elements and a convective turnover time of
$\sim$100-300 days, we expect to wait $\sim$30-100 years between major mass
ejection events.

Each mass loss event would last for a convective turnover time of $\sim$100
days, ejecting mass at a rate $\sim\\!10^{-4}-10^{-3}\,M_{\odot}$/yr (see
green line in Figure 4), ejecting a total mass of $M_{\rm ej}\\!\sim\\!3\times
10^{-5}-3\times 10^{-4}\,M_{\odot}$ and about $3\times 10^{-7}-3\times
10^{-6}\,M_{\odot}$ of dust. For a dust-laden gas opacity of $\kappa_{\rm
d}\sim 3\,{\rm cm}^{2}/{\rm g}$, this mass can have an optical depth $\tau\sim
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_{\odot}/{\rm 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_{\rm d}/v_{\rm
wind}\sim 2\,{\rm yr}$, the probability of observing such a star with a very
dense CSM at the moment of a SN explosion is near unity.

###  5.5 Comparison with supernovae and SN 2023ixf

Observations of type-IIP supernovae have frequently indicated the presence of
dense CSM within $\sim\\!10^{15}$ cm ($\sim$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 $\sim$5 days after shock-breakout (Bruch et al., 2023) as this CSM is
swept up by the SN ejecta that travels outwards at $\sim\\!10^{9}$ cm/s
($\sim\\!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 $\sim$10
stellar radii ($\sim 5\times 10^{14}\,{\rm cm}$) of exploding RSGs. This
material will be swept up by the SN ejecta in $\sim$5 days, in line with the
disappearance of flash-ionized 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 $\dot{M}\sim
10^{-5}M_{\odot}/{\rm 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 $\sim\\!10^{-14}\,{\rm g}/{\rm
cm}^{3}$ at $\sim$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 $\sim$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,
$\rho\sim\\!5\times 10^{-14}\,{\rm g}/{\rm cm}^{3}$, is slightly larger than
our time-averaged model, which has a volume-averaged CSM density of
$\rho\sim\\!5\times 10^{-15}\,{\rm g}/{\rm 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 $\sim\\!2\times 10^{14}\,{\rm 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
$\dot{M}\sim 10^{-4}\,M_{\odot}/{\rm 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 $\sim\\!10^{15}\,{\rm cm}$ is inferred to have a density of
$\sim\\!4\times 10^{-16}\,{\rm g}/{\rm 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 $\sim\\!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 $\rho\sim 10^{-12}\,{\rm cm}$ at $\sim$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 $\sim\\!10^{-2}-10^{-1}\,M_{\odot}/{\rm 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 V-band 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 & 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), in contrast to most RSGs and SN progenitors. Such pulsation is thought
to be crucial to the mass loss of AGB stars by driving shocks through the
atmosphere (Höfner & Olofsson, 2018; Freytag et al., 2017). Our model could
possibly account for this by replacing the convective velocity $v_{\rm con}$
with the pulsational velocity near the photosphere, which can be supersonic
and hence larger than $v_{\rm 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 & Cantiello, 2010).

###  5.6 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\gtrsim 2v_{\rm 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 &
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 log-normal 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 $\kappa\sim
10^{-4}-10^{-2}\,{\rm cm^{2}}/{\rm g}$ (Plez et al., 1992). The corresponding
Eddington factor, i.e., the ratio between radiative acceleration and gravity,
is $f_{\rm ed}=g/g_{\rm rad}=L/L_{\rm ed}=L\kappa/(4\pi GMc)\sim
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 &
Plez, 2007). Since our predicted mass loss rates are strongly dependent on
$v_{\rm 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 & Heyl (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 & 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.

##  6 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 $\sim$5-10 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_{\rm con}^{2}/g\ll R$, which would cause a steeply decreasing density for
shock speeds comparable to the surface convective speeds $v_{\rm 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 3.1). 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_{\rm con}/v_{\rm 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 $\sim
10^{-6}\,M_{\odot}/{\rm yr}$ for $15-20\,M_{\odot}$ RSGs, ii) substantially
lower mass loss rates for less massive RSGs as found by recent observations,
iii) the absence of RSGs for luminosities $L\gtrsim 3\times 10^{5}\,L_{\odot}$
due to rapid boil-off of their hydrogen envelopes. Our models also predict
pre-supernova mass loss rates of $\sim\\!4\times 10^{-7}-4\times
10^{-5}\,M_{\odot}$/yr, a factor 10-100 times larger than mass loss rates
during core helium burning.

Our models predict chromospheric masses of $\sim\\!10^{-2}\,M_{\odot}$, 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 constant-velocity 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_{\rm 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.

##  7 Acknowledments

We thank Luc Dessart and Jared Goldberg for useful discussions, and Kostas
Antoniadis for supplying mass loss measurements. This work benefited from
collaborative meetings supported by the Gordon and Betty Moore Foundation
(grant GBMF5076). This research was supported by the Munich Institute for
Astro-, Particle and BioPhysics (MIAPbP) which is funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany´s
Excellence Strategy – EXC-2094 – 390783311. This research is also supported by
the Sherman Fairchild Postdoctoral Fellowship at the California Institute of
Technology.

## References

  * Antoniadis et al. (2024) Antoniadis K., Bonanos A. Z., de Wit S., Zapartas E., Munoz-Sanchez G., Maravelias G., 2024, arXiv e-prints, p. arXiv:2401.15163
  * Arroyo-Torres et al. (2015) Arroyo-Torres B., et al., 2015, A&A, 575, A50
  * Beasor & Smith (2022) Beasor E. R., Smith N., 2022, ApJ, 933, 41
  * Beasor et al. (2020) Beasor E. R., Davies B., Smith N., van Loon J. T., Gehrz R. D., Figer D. F., 2020, MNRAS, 492, 5994
  * Bertschinger & Chevalier (1985) Bertschinger E., Chevalier R. A., 1985, ApJ, 299, 167
  * Bostroem et al. (2023) Bostroem K. A., et al., 2023, ApJ, 956, L5
  * Bowen (1988) Bowen G. H., 1988, ApJ, 329, 299
  * Bruch et al. (2021) Bruch R. J., et al., 2021, ApJ, 912, 46
  * Bruch et al. (2023) Bruch R. J., et al., 2023, ApJ, 952, 119
  * Chandrasekhar (1934) Chandrasekhar S., 1934, MNRAS, 94, 444
  * Chiavassa et al. (2011a) Chiavassa A., et al., 2011a, A&A, 528, A120
  * Chiavassa et al. (2011b) Chiavassa A., Freytag B., Masseron T., Plez B., 2011b, A&A, 535, A22
  * Chiavassa et al. (2024) Chiavassa A., Kravchenko K., Goldberg J. A., 2024, arXiv e-prints, p. arXiv:2402.00187
  * Coughlin et al. (2018) Coughlin E. R., Quataert E., Ro S., 2018, ApJ, 863, 158
  * Cranmer & Saar (2011) Cranmer S. R., Saar S. H., 2011, ApJ, 741, 54
  * Davies et al. (2018) Davies B., Crowther P. A., Beasor E. R., 2018, MNRAS, 478, 3138
  * Davies et al. (2022) Davies B., Plez B., Petrault M., 2022, MNRAS, 517, 1483
  * Dessart & Jacobson-Galán (2023) Dessart L., Jacobson-Galán W. V., 2023, A&A, 677, A105
  * Dessart et al. (2017) Dessart L., Hillier D. J., Audit E., 2017, preprint,  (arXiv:1704.01697)
  * Dupree et al. (2022) Dupree A. K., et al., 2022, ApJ, 936, 18
  * Field (1974) Field G. B., 1974, ApJ, 187, 453
  * Freytag & Höfner (2008) Freytag B., Höfner S., 2008, A&A, 483, 571
  * Freytag & Höfner (2023) Freytag B., Höfner S., 2023, A&A, 669, A155
  * Freytag et al. (2017) Freytag B., Liljegren S., Höfner S., 2017, A&A, 600, A137
  * Fuller (2017) Fuller J., 2017, MNRAS, 470, 1642
  * Goldberg et al. (2022a) Goldberg J. A., Jiang Y.-F., Bildsten L., 2022a, ApJ, 929, 156
  * Goldberg et al. (2022b) Goldberg J. A., Jiang Y.-F., Bildsten L., 2022b, ApJ, 933, 164
  * Gräfener & Vink (2016) Gräfener G., Vink J. S., 2016, MNRAS, 455, 112
  * Grefenstette et al. (2023) Grefenstette B. W., Brightman M., Earnshaw H. P., Harrison F. A., Margutti R., 2023, The Astrophysical Journal Letters, 952, L3
  * Haynie & Piro (2021) Haynie A., Piro A. L., 2021, ApJ, 910, 128
  * Hill & Willson (1979) Hill S. J., Willson L. A., 1979, ApJ, 229, 1029
  * Hiramatsu et al. (2021) Hiramatsu D., et al., 2021, Nature Astronomy, 5, 903
  * Hiramatsu et al. (2023) Hiramatsu D., et al., 2023, ApJ, 955, L8
  * Höfner & Olofsson (2018) Höfner S., Olofsson H., 2018, A&A Rev., 26, 1
  * Höfner et al. (2003) Höfner S., Gautschy-Loidl R., Aringer B., Jørgensen U. G., 2003, A&A, 399, 589
  * Hopkins (2013) Hopkins P. F., 2013, MNRAS, 430, 1880
  * Hosseinzadeh et al. (2023) Hosseinzadeh G., et al., 2023, ApJ, 953, L16
  * Hoyle & Wickramasinghe (1962) Hoyle F., Wickramasinghe N. C., 1962, MNRAS, 124, 417
  * Irani et al. (2023) Irani I., et al., 2023, arXiv e-prints, p. arXiv:2310.16885
  * Ireland et al. (2008) Ireland M. J., Scholz M., Wood P. R., 2008, MNRAS, 391, 1994
  * Jacobson-Galán et al. (2022) Jacobson-Galán W. V., et al., 2022, ApJ, 924, 15
  * Jacobson-Galán et al. (2023) Jacobson-Galán W. V., et al., 2023, ApJ, 954, L42
  * Jadlovský et al. (2023) Jadlovský D., Krtička J., Paunzen E., Štefl V., 2023, New A, 99, 101962
  * Jencson et al. (2023) Jencson J. E., et al., 2023, ApJ, 952, L30
  * Josselin & Plez (2007) Josselin E., Plez B., 2007, A&A, 469, 671
  * Kee et al. (2021) Kee N. D., Sundqvist J. O., Decin L., de Koter A., Sana H., 2021, A&A, 646, A180
  * Khazov et al. (2016) Khazov D., et al., 2016, ApJ, 818, 3
  * Khouri et al. (2024) Khouri T., et al., 2024, arXiv e-prints, p. arXiv:2402.13676
  * Kilpatrick et al. (2023) Kilpatrick C. D., et al., 2023, ApJ, 952, L23
  * Levesque et al. (2005) Levesque E. M., Massey P., Olsen K. A. G., Plez B., Josselin E., Maeder A., Meynet G., 2005, ApJ, 628, 973
  * Li et al. (2023) Li G., et al., 2023, arXiv e-prints, p. arXiv:2311.14409
  * López Ariste et al. (2022) López Ariste A., Georgiev S., Mathias P., Lèbre A., Wavasseur M., Josselin E., Konstantinova-Antova R., Roudier T., 2022, A&A, 661, A91
  * Ma et al. (2023) Ma J.-Z., Chiavassa A., de Mink S. E., Valli R., Justham S., Freytag B., 2023, arXiv e-prints, p. arXiv:2311.16885
  * Marshall et al. (2004) Marshall J. R., van Loon J. T., Matsuura M., Wood P. R., Zijlstra A. A., Whitelock P. A., 2004, MNRAS, 355, 1348
  * Matzner & McKee (1999) Matzner C. D., McKee C. F., 1999, ApJ, 510, 379
  * Mauron & Josselin (2011) Mauron N., Josselin E., 2011, A&A, 526, A156
  * Meza Retamal et al. (2024) Meza Retamal N., et al., 2024, arXiv e-prints, p. arXiv:2401.04027
  * Montargès et al. (2021) Montargès M., et al., 2021, Nature, 594, 365
  * Moriya et al. (2017) Moriya T. J., Yoon S.-C., Gräfener G., Blinnikov S. I., 2017, preprint,  (arXiv:1703.03084)
  * Moriya et al. (2018) Moriya T. J., Förster F., Yoon S.-C., Gräfener G., Blinnikov S. I., 2018, MNRAS, 476, 2840
  * Morozova et al. (2015) Morozova V., Piro A. L., Renzo M., Ott C. D., Clausen D., Couch S. M., Ellis J., Roberts L. F., 2015, ApJ, 814, 63
  * Morozova et al. (2017) Morozova V., Piro A. L., Valenti S., 2017, ApJ, 838, 28
  * Morozova et al. (2018) Morozova V., Piro A. L., Valenti S., 2018, ApJ, 858, 15
  * Morozova et al. (2020) Morozova V., Piro A. L., Fuller J., Van Dyk S. D., 2020, ApJ, 891, L32
  * Moseley et al. (2019) Moseley E. R., Squire J., Hopkins P. F., 2019, MNRAS, 489, 325
  * Nakar & Sari (2010) Nakar E., Sari R., 2010, ApJ, 725, 904
  * Ofek et al. (2014) Ofek E. O., et al., 2014, ApJ, 789, 104
  * Ohnaka et al. (2017) Ohnaka K., Weigelt G., Hofmann K. H., 2017, Nature, 548, 310
  * Paxton et al. (2011) Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton B., et al., 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton B., et al., 2015, ApJS, 220, 15
  * Plez et al. (1992) Plez B., Brett J. M., Nordlund A., 1992, A&A, 256, 551
  * Rabatin & Collins (2023a) Rabatin B., Collins D. C., 2023a, MNRAS, 521, L64
  * Rabatin & Collins (2023b) Rabatin B., Collins D. C., 2023b, MNRAS, 525, 297
  * Reimers (1975) Reimers D., 1975, Memoires of the Societe Royale des Sciences de Liege, 8, 369
  * Sakurai (1960) Sakurai A., 1960, Communs. Pure and Appl. Math., 13
  * Schirrmacher et al. (2003) Schirrmacher V., Woitke P., Sedlmayr E., 2003, A&A, 404, 267
  * Schröder & Cuntz (2005) Schröder K. P., Cuntz M., 2005, ApJ, 630, L73
  * Shu (1992) Shu F. H., 1992, The physics of astrophysics. Volume II: Gas dynamics.
  * Sjouwerman et al. (1998) Sjouwerman L. O., van Langevelde H. J., Winnberg A., Habing H. J., 1998, A&AS, 128, 35
  * Smith et al. (2023) Smith N., Pearson J., Sand D. J., Ilyin I., Bostroem K. A., Hosseinzadeh G., Shrestha M., 2023, ApJ, 956, 46
  * Soker (2021) Soker N., 2021, ApJ, 906, 1
  * Squire & Hopkins (2018) Squire J., Hopkins P. F., 2018, MNRAS, 477, 5011
  * Steinwandel et al. (2022) Steinwandel U. P., Kaurov A. A., Hopkins P. F., Squire J., 2022, MNRAS, 515, 4797
  * Strotjohann et al. (2021) Strotjohann N. L., et al., 2021, ApJ, 907, 99
  * Svirski et al. (2012) Svirski G., Nakar E., Sari R., 2012, ApJ, 759, 108
  * Taniguchi et al. (2022) Taniguchi D., Yamazaki K., Uno S., 2022, Nature Astronomy, 6, 930
  * Thirumalai & Heyl (2010) Thirumalai A., Heyl J. S., 2010, MNRAS, 409, 1669
  * Thirumalai & Heyl (2012) Thirumalai A., Heyl J. S., 2012, MNRAS, 422, 1272
  * Tsuna et al. (2021) Tsuna D., Kashiyama K., Shigeyama T., 2021, ApJ, 914, 64
  * Tsuna et al. (2023a) Tsuna D., Murase K., Moriya T. J., 2023a, ApJ, 952, 115
  * Tsuna et al. (2023b) Tsuna D., Murase K., Moriya T. J., 2023b, ApJ, 952, 115
  * Vasylyev et al. (2023) Vasylyev S. S., et al., 2023, ApJ, 955, L37
  * Vink & Sabhahit (2023) Vink J. S., Sabhahit G. N., 2023, A&A, 678, L3
  * Wedemeyer et al. (2017) Wedemeyer S., Kučinskas A., Klevas J., Ludwig H.-G., 2017, A&A, 606, A26
  * Wen et al. (2024) Wen J., Gao J., Yang M., Chen B., Ren Y., Wang T., Jiang B., 2024, AJ, 167, 51
  * Whitham (2011) Whitham G., 2011, Linear and Nonlinear Waves. Pure and Applied Mathematics: A Wiley Series of Texts, Monographs and Tracts, Wiley, https://books.google.com/books?id=84Pulkf-Oa8C
  * Winters et al. (1997) Winters J. M., Fleischer A. J., Le Bertre T., Sedlmayr E., 1997, A&A, 326, 305
  * Woosley & Heger (2015) Woosley S. E., Heger A., 2015, ApJ, 810, 34
  * Yang et al. (2023) Yang M., et al., 2023, A&A, 676, A84
  * Yaron et al. (2017) Yaron O., et al., 2017, preprint,  (arXiv:1701.02596)
  * Yoon & Cantiello (2010) Yoon S.-C., Cantiello M., 2010, ApJ, 717, L62
  * Zhang et al. (2023) Zhang J., et al., 2023, arXiv e-prints, p. arXiv:2309.01998
  * Zimmerman et al. (2023) Zimmerman E. A., et al., 2023, arXiv e-prints, p. arXiv:2310.10727
  * de Jager et al. (1988) de Jager C., Nieuwenhuijzen H., van der Hucht K. A., 1988, A&AS, 72, 259
  * van Loon et al. (2001) van Loon J. T., Zijlstra A. A., Bujarrabal V., Nyman L. Å., 2001, A&A, 368, 950
  * van Loon et al. (2005) van Loon J. T., Cioni M. R. L., Zijlstra A. A., Loup C., 2005, A&A, 438, 273