∎

11institutetext: A. Chiavassa 22institutetext: Université Côte d’Azur,
Observatoire de la Côte d’Azur, CNRS, Lagrange, CS 34229, Nice, France
and
Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, 85741
Garching, Germany
22email<EMAIL_ADDRESS>33institutetext: K. Kravchenko
44institutetext: Max Planck Institute for Extraterrestrial Physics,
Gießenbachstraße 1, 85748 Garching, Germany 55institutetext: J. A. Goldberg
66institutetext: Center for Computational Astrophysics, Flatiron Institute,
162 5th Avenue, 10010, New York, NY, USA

# Signatures of convection in the atmospheres of cool evolved stars

Andrea Chiavassa  Kateryna Kravchenko  Jared A. Goldberg

(Received: date / Accepted: date)

###### Abstract

Evolved cool stars of various masses are major cosmic engines, delivering
substantial mechanical and radiative feedback to the interstellar medium
through strong stellar winds and supernova ejecta. These stars play a pivotal
role in enriching the interstellar medium with vital chemical elements that
constitute the essential building blocks for the formation of subsequent
generations of stars, planets, and potentially even life.
Within the complex tapestry of processes occurring in the atmospheres of these
cool and luminous stars, convection takes center stage. Convection is a non-
local, complex phenomenon marked by non-linear interactions across diverse
length scales within a multi-dimensional framework. For these particular
stars, characterized by their considerable luminosities and extensive scale
heights, convection transitions to a global scale. This transition is
facilitated by the transmission of radiative energy through the non-uniform
outer layers of their atmospheres.
To have a full understanding of this phenomenon, the application of global
comprehensive 3D radiation-hydrodynamics simulations of stellar convection is
of paramount importance. We present two state-of-the-art numerical codes:
CO5BOLD and Athena++. Furthermore, we provide a view on their applications as:
pivotal roles in enabling a comprehensive investigation into the dynamic
processes linked to convection; and critical tools for accurately modeling the
emissions produced during shock breakouts in Type II-P Supernovae.

###### Keywords:

Stars: atmospheres hydrodynamics radiative transfer methods: numerical

###### pacs:

85A30 85A25

##  1 Introduction

Evolved cool stars of various masses are major cosmic engines, providing
strong mechanical and radiative feedback on their host environment (Langer,
2012). Through strong stellar winds and supernova ejecta, they enrich the
interstellar medium with chemical elements, which are the building blocks of
planets and life. In particular, these objects are known to propel strong
stellar winds at speeds varying with stellar brightness, evolutionary phase,
and chemical composition. Yet, stellar evolution models are not able to
reproduce these winds without ”ad hoc” physics. Therefore, a full
understanding of stellar evolution in the near and distant Universe, and its
impact on the cosmic environment, cannot be achieved without a detailed
knowledge of wind physics across the life cycle of these stars. This requires
tracing the total mass ejected as well as its composition, the velocity of the
winds, and the overall geometry of the circumstellar envelope. Evolved cool
stars are objects that have reached the late phases of their evolution when
the nuclear fuel in the interior is almost used up. These stars grow
dramatically in size and become Asymptotic Giant Branch (AGB, the initial mass
of $0.8<M<8\,M_{\odot}$, the exact value of the upper limit depends on the
treatment of convection) or Red Supergiant (RSG, $M>8\,M_{\odot}$) stars. AGBs
are characterized by a very low effective temperature (lower than 3000 K) and
a high luminosity (about 100–1000 $L_{\odot}$) which correspond to stellar
radii of several hundred solar radii. Their high rate of mass-loss rate
($\dot{M}=10^{-6}$–$10^{-4}\,M_{\odot}/\mathrm{yr}$, Wood et al., 1992) are
the results of a wind acceleration mechanism based on radiation pressure on
microscopic solid-state particles (dust grains), formed in atmospheres that
are levitated by pulsation. State-of-the-art models are still insufficient to
provide a full and consistent theory of their mass-loss (Höfner and Olofsson,
2018). RSGs are precursors of core-collapse supernovae with high luminosities
($L>1000\,L_{\odot}$) with effective temperatures between 3450 and 4100 K and
stellar radii up to several hundreds of $R_{\odot}$, or even more than
$1000\,R_{\odot}$ (Levesque et al., 2005). Several mechanisms triggering mass-
loss have been discussed, including magneto-hydrodynamic waves (e.g. Cranmer
and Saar 2011), turbulent pressure from trans-sonic convection (e.g. Kee et
al. 2021), and radiation pressure on molecules and dust (Josselin and Plez,
2007), but still there is no realistic quantitative wind model (Meynet et al.,
2015) that can explain the observed broad range of mass-loss rates
($\dot{M}=10^{-7}$–$10^{-4}\,M_{\odot}/\mathrm{yr}$, De Beck et al., 2010).

A full picture of all the physical processes that simultaneously trigger and
shape the convective envelopes and strong winds of evolved cool stars is still
missing. Indeed, the mass-loss mechanism is hard to discern because it
involves a range of interacting, time-dependent physical processes on
microscopic and macroscopic scales coupled with dynamical phenomenon such as
convection and pulsation in sub-photospheric layers, strong radiating shocks
in the atmosphere, and dust condensation as well as radiative acceleration in
the wind forming regions (Höfner and Olofsson, 2018). By depositing angular
momentum, (sub)stellar companions are known to shape the outflow of cool
evolved stars (Decin et al., 2020).

In this framework, three physical ingredients should play a predominant role,
depending on the stellar type, in initiating and maintaining the high mass
loss. In the first place, the evolution of these objects is impacted by
stellar convection. Convection is a difficult process to understand because it
is non-local, three-dimensional, and involves non-linear interactions over
many disparate length scales, and it is partially responsible for transporting
heat up to the visible surface (Nordlund et al., 2009). In AGB and RSG
atmospheres, convection is inferred from a few giant structures observed at
the stellar surface with sizes comparable to the stellar radius and evolving
on weekly or yearly time scales (e.g. Ma et al. 2023; Montargès et al. 2021;
Paladini et al. 2018a; Montargès et al. 2018; Chiavassa et al. 2010, 2011b).
Such large-scale convective structures are expected theoretically due to the
convective plume size approximately scaling with the larger pressure scale-
heights in these giant envelopes (e.g. Böhm-Vitense, 1958; Trampedach et al.,
2014; Chiavassa et al., 2011a; Freytag et al., 2017; Goldberg et al., 2022a).
This results in more extreme atmospheric conditions than in the Sun: very
large variations in velocity, density, and temperature produce strong
radiative shocks in their extended atmosphere that can cause the gas to
levitate and thus contribute to mass-loss (Chiavassa et al., 2011a; Freytag et
al., 2017; Höfner and Freytag, 2019; Ahmad et al., 2023). In the outer layers
of these extreme atmospheres, radiation begins to carry non-negligible amounts
of the flux, yet the the opacity, set by the (varied) density and temperature,
can experience order-of-magnitude fluctuations from H and He opacity peaks.
Because the luminosities are so high, these fluctuations in opacity will
affect the dynamics of the convective envelope.

The second ingredient is dust. In the case of the slow and massive winds of
AGB stars, the wind is commonly assumed to be driven by radiative pressure on
dust grains, forming in the extended stellar atmospheres Höfner and Olofsson
(2018). The pulsations, together with large-scale convective flows, trigger
strong radiative shock waves in the stellar atmospheres, which intermittently
lift gas to distances where temperatures are low enough to let dust form
Höfner and Freytag (2019). This has been recently shown in the theoretical
work of Freytag and Höfner (2023), where the simulations are capable to follow
the flow of matter from the stellar interior, through the dynamical
atmosphere, and into the wind-acceleration region, while taking non-
equilibrium dust formation and the interaction of matter with radiation into
account. The RSG phase may also be characterised by the presence of dust close
to a star (Haubois et al., 2019; Montargès et al., 2021) as well as the
inferences of circumstellar matter in the pulsating progenitor of a Type II
supernova (Kozyreva et al., 2022). A similar picture as for AGB stars could
also be true for RSGs where pulsation-driven winds may dominate (e.g. Yoon and
Cantiello 2010).

The third ingredient is the magnetic field. Cranmer and Saar (2011) presented
a predictive description of mass-loss, based on Alfvén-wave-driven wind.
Alfvén waves (i.e., incompressible, transverse magneto-hydrodynamic waves) are
expected to provide a source of heating by dissipation. It should be noted
that Alfvén-wave-driven winds require open field lines, radially directed away
from the star, in order for the gas to be accelerated and escape (Höfner and
Olofsson, 2018). RSGs and AGBs have low-intensity integrated magnetic fields
(of the order of 1-10 Gauss) that have been identified and monitored over
several years (Tessore et al., 2017; Lèbre et al., 2014; Aurière et al.,
2010). However, the knowledge of the magnetic field strength and geometry is
still limited and its origin from a possible astrophysical dynamo in these
stars would most likely be very different from the dynamo at work in solar-
type stars due to the small radiative core, slow rotation and the fact that
only a few convection cells are present at their surface at any given time
(Aurière et al., 2010; Freytag et al., 2002).

To provide quantitative constraints on the physics of cool evolved stars,
observational techniques have reached such a level of excellence that it is
now possible to reconstruct spatially resolved images of stellar surfaces in
the near IR with interferometric (e.g. Chiavassa et al. 2022a, 2020; Montargès
et al. 2018; Paladini et al. 2018a; Ohnaka et al. 2017; Haubois et al. 2009)
or astrometric techniques (Chiavassa et al., 2022b, 2018b). Simultaneous
intensity and polarization stellar surface images have been obtained at
unprecedented resolutions with SPHERE@VLT (e.g. Kervella et al. 2016). Images
obtained with ALMA were used to measure the rotation of the molecular envelope
of the RSG star Betelgeuse (Kervella et al., 2018). Finally,
spectropolarimetric observations (with the instruments Narval-TBL and
ESPaDOnS-CFHT) have made it possible, among other things, to detect and
characterize the surface magnetism and dynamics of RSGs (López Ariste et al.,
2022; Mathias et al., 2018; Aurière et al., 2010) and AGBs (López Ariste et
al., 2019; Lèbre et al., 2014).

In the following Sections, we concentrate on the stellar atmospheric dynamics
based on convection-related processes, leaving further reviews on dust- or
magneto-related subjects for later works.

##  2 3D radiation-hydrodynamics simulations of stellar convection

The interpretation of high-resolution spatially and spectrally resolved
observations of evolved stars requires realistic modeling that takes into
account most of the processes at work in the atmosphere (e.g., convection,
shocks, pulsation, radiative transfer, ionization, molecules and dust
formation, magnetic fields).

Despite the apparent sphericity of stars, all the physical processes happening
in their atmospheres are truly three-dimensional and therefore, they require a
three-dimensional treatment, unlike the traditional one-dimensional approach
in stellar-evolution modeling. For instance, one crucial point is the
radiative transfer treatment that has to be simultaneously solved through a
highly turbulent medium with large density variations over optical depths
ranging from far above unity down to the photospheric regions where the
radiation’s contribution is larger. Another related point is the convection-
related processes, that involve non-linear interactions over many disparate
length scales. In this context, the use of numerical 3D radiation-
hydrodynamics simulations of stellar convection is extremely important. In
recent years, with increased computational power, it has been possible to
compute grids of 3D radiation-hydrodynamics (RHD) simulations of the whole
stellar envelope that are used to predict reliable synthetic spectra and
images for several stellar types. The first realistic simulations of solar
granulation were performed by Nordlund (1982) and included 3D, time-dependent
hydrodynamics (with moderate spatial resolution) and non-local radiative
energy transfer, already then with a simple treatment of the frequency-
dependence of the opacities. After this pioneering work and more recently
several 3D simulations (in ”box-in-a-star” configuration) grids have been
produced for a substantial portion of the Hertzsprung-Russell (Magic et al.,
2013; Trampedach et al., 2013; Ludwig et al., 2009). They are used to obtain
precise determinations of stellar parameters (e.g. Bigot et al. 2011; Creevey
et al. 2012), radial velocity (e.g. Allende Prieto et al. 2013; Chiavassa et
al. 2018a), chemical abundance (e.g. Asplund et al. 2009; Caffau et al. 2011;
Magg et al. 2022), photometric colours (e.g. Bonifacio et al. 2018; Chiavassa
et al. 2018a), and on planet detection/characterisation (e.g. Magic et al.
2015a; Chiavassa and Brogi 2019; Cegla 2019). Moreover, these grids have
yielded a physical understanding of the nature of RHD convection, proving a
quantitative way to set the outer boundary condition in 1D stellar models
(Trampedach et al., 2014; Salaris and Cassisi, 2015; Magic and Weiss, 2016;
Mosumgaard et al., 2018; Sonoi et al., 2019).

However, the large sizes of the turbulent, extended outer envelope of evolved
stars (such as AGB and RSG) manifest inherently 3D convective properties that
need the use of global (i.e., including a major part of the convective
envelope as well as the near environment of the star) and CPU-intensive
simulations. In these stars, the convection-related structures are comparable
to the stellar diameter and need a large enough domain. To account for that,
global ”star-in-a-box” simulations are used: they include the majority of the
convective envelope as well as the near environment of the star and differ
essentially in boundary conditions and the gravitational potential from local
”box-in-a-star” models that use constant surface gravity (i.e. plane-parallel
approximation) and simulate only a relatively small area and shallow depths of
a star with parameters close to the main sequence and Red Giant Branch stars.
The ”star-in-a-box” simulations return very large surface contrast, due to the
combination of violent convective flows and high-amplitude waves, which is
very different from the solar case of both smaller convective structures and
lower mode amplitudes (Freytag et al., 2012). There are two numerical codes
that are able to perform such simulations: (i) the CO5BOLD code (Freytag et
al., 2012) has been used to simulate AGB (Freytag et al., 2017; Ahmad et al.,
2023) and RSG (Chiavassa et al., 2011a) stars; (ii) the Athena++ code (Stone
et al., 2020) that can simulate RSG stars (Goldberg et al., 2022a) in addition
to other luminous massive stars, (Jiang et al., 2018; Schultz et al., 2022,
2023). With these two numerical approaches, it is possible to follow the
flowing matter in complete 3D geometry: from the deep internal layers, where
non-stationary convection produces sound waves; to the surface, where the
otherwise standing waves can partially escape from the acoustic cavity of the
interior, and speed up due to the exponential decrease in density and turn
into shocks, that propagates all the way to the outer circumstellar envelope.

###  2.1 CO5BOLD and Athena++ codes

CO5BOLD solves the coupled equations of compressible hydrodynamics and non-
local radiative energy transport in the presence of a fixed external
spherically symmetric gravitational field on a 3D cartesian grid. It is
important to include an accurate gravitational potential because the mass in
the simulation domain is not negligible and the envelope is loosely bound. The
hydrodynamics scheme of CO5BOLD uses a finite-volume approach. By means of
operator (directional) splitting (Strang, 1968), the 2D or 3D problem is
reduced to one dimension. To compute the fluxes across each cell boundary in
every 1D column direction, an approximate 1D Riemann solver of Roe type (Roe,
1986) is applied, including appropriate modifications to account an adapted
Equation of State, a non-equidistant grid, and a gravitational potential
(Freytag et al., 2012). CO5BOLD can also handle the calculation of the
magneto-hydrodynamic equations (Riva and Steiner, 2022) using a flux-
conserving HLL solver (Harten et al., 1983). Finally, CO5BOLD may account for
rotation performing calculations in a corotating frame (Freytag et al., 2017):
a centrifugal potential is added to the gravitational potential resulting in
an application of the Coriolis force at each hydrodynamics step. All six
surfaces of the computational box employ the same open boundary condition. The
radiative transport is performed using a short-characteristics method (Freytag
et al., 2012) using Rosseland mean opacity tables (merged at around 12 000 K
from high-temperature OPAL data and low-temperature PHOENIX data, Iglesias et
al., 1992; Hauschildt et al., 1997, respectively). Different works report
global simulations employing three or five opacity bins (Chiavassa et al.,
2011a; Höfner and Freytag, 2019; Freytag and Höfner, 2023).

Athena++ is an adaptive mesh refinement code solving the hydrodynamic
equations based on a task-based dynamic execution (similar to what has been
implemented in other codes such as DISPATCH, Nordlund et al., 2018). To
provide robust and accurate shock capturing, a Godunov-type method is used to
ensure the non-oscillatory spatial reconstruction of the fluid variables at
cell interfaces, a Riemann solver to compute upwind fluxes and electric fields
at cell faces, and a time-integration algorithm to advance the solution (Stone
et al., 2020). Athena++ uses distributed-memory parallelism through domain
decomposition and is capable of employing both non-relativistic and
relativistic magnetohydrodynamics. The simulations of RGS stars use a
spherical polar coordinate system with 256 bins in azimuthal angle and 128
uniform bins in polar angle, assuming a non-rotating stellar model and
neglecting magnetic fields (Goldberg et al., 2022a). The code includes an
accurate spherically symmetric gravitational potential $GM(r)/r$ which is a
function of both the mass inside the simulation domain’s inner boundary as
well as the mass in the simulated envelope ($M(r)=M_{\rm IB}+\Delta m(r)$
where $\Delta m$ is the total enclosed mass at radius $r$ from summing each
successive shell). The radiative transport is fully coupled to the
hydrodynamics in Athena++ directly solving the time-dependent, frequency-
integrated (grey) radiative transport equation for specific intensities over
discrete ordinates (Jiang, 2021) with 120 angles. Radiation is solved for in
the co-moving frame and then self-consistently accounted for as a source term
in both the energy and momentum equations by transforming to the lab frame.
The Rosseland and Planck mean absorption opacities are determined by
interpolation of the OPAL opacity tables (Iglesias and Rogers, 1996). Similar
numerical setups have been successfully used to model both local patches and
global simulations of stellar envelopes in different locations of the HR
diagram (Jiang et al., 2015, 2018; Schultz et al., 2022, 2023). The ability to
handle multi-group opacities exists in Athena++ (Jiang, 2022), but is
presently untested for stellar envelopes.

###  2.2 Brief summary of results and mutual comparison

The first CO5BOLD simulation were published by Freytag et al. (2002) and,
since then, have proved to be realistic and indispensable tools for the
understanding of the various complex dynamical processes in evolved stars.
Since 2010, CO5BOLD simulations have been used to interpret the large majority
of interferometric imaging observations (for a short summary of imaging
programs on stellar convection, see Paladini et al., 2018b) with the VLTI and
CHARA telescopes: in fact, large convection-related structures and large
intensity contrasts render these features observable in detail with current
interferometers. More recently, CO5BOLD simulations have also been used to
interpret convection cycles using spectroscopy (e.g. Kravchenko et al. 2021
and Section 4.1).
While CO5BOLD simulations have been available since Freytag et al. (2002), the
employment of Athena++ to compute RSG stars is very recent and and already
important results have been reported. Goldberg et al. (2022a) showed an
unusual outcome of an inverse correlation between density and radial velocity
in the outer layers, uncharacteristic of typical stellar convection, and a
density structure much more extended than 1D stellar evolution predictions.
The authors also managed to calibrate the mixing length ($\alpha$ = 3 – 4) for
the mixing length theory (MLT), showing that the radiation pressure provides
$\approx$ 1/3 of the support against gravity in this region. The extended and
inhomogeneous atmosphere has also an impact on the shock breakout of Type II-P
supernovae, as shown by Goldberg et al. (2022b): the large-scale order-of-
magnitude 3D fluctuations in density cause the shock to break out at different
radii at different times, and the duration of the shock breakout is prolonged
over $\approx$3 – 6 hr, smeared out by the shock-crossing time of the
inhomegeneous outer layers. The authors conclude that the intrinsically 3D
nature of the envelope precludes the possibility of using observed UV rise
times to measure the stellar radius via light-travel time effects, but enables
shock-breakout observations to probe the scale of the convective
inhomogeneities when combined with additional information about the shock
velocity from other supernova observables.
A quantitative comparison between the results of the two codes is reported in
Goldberg et al. (2022a), with a few points: (i) CO5BOLD code neglects the
radiation pressure in deeper layers of the star (an approximation that is only
acceptable for low optical depth, where radiative transport dominates and
where the observed flux is formed), inhibiting the ability to correctly
simulate the deeper nearly-constant-entropy convective zone there. As a
consequence, the CO5BOLD RSG simulations exhibit a positive entropy gradient
at about $\approx$75$\%$ of the stellar radius (Chiavassa et al., 2011a) while
Athena++ simulation predicts a nearly-flat, declining entropy profile
($ds/dr\lesssim 0$) throughout the convective envelope more consistent with 1D
models of stellar convection. (ii) Both codes show coherent convective plumes
across the star with comparable radial velocities of tens of km/s and
comparable convective overturn times of hundreds of days. (iii) Simulations
with both codes display high-entropy material above the expected stellar
photospheric radii with density fluctuations developing far from the stellar
photosphere. (iv) Both codes predict that convection contributes significantly
to the energy transport in the stellar interior and (v) that the turbulent
energy density from the vigorous convective motions dominates over the thermal
energy in the region above the photospheric stellar radius.

##  3 The dynamic atmosphere of evolved stars as seen by state-of-the-art
numerical simulations

###  3.1 The dynamical atmosphere as seen by CO5BOLD

In this section, we concentrate on the results obtained with CO5BOLD to report
the outcome of the convective cycles extracted from high spectral resolution
spectroscopy (Section 4.1). We present in detail the dynamical processes
emerging from the simulations of AGB and RSG stellar atmospheres (Table 1).

####  3.1.1 Starting a simulation

As described in Chiavassa et al. (2011a), the principal requirements that
determine the type of the simulated star are (i) the input luminosity into the
core, (ii) the stellar mass (i.e., the gravitational potential), (iii) the
chemical abundances that affect the atomic physics through the opacity and
equation of state tables.
Moreover, several additional parameters influence the outcome of a simulation
to some extent: the formulation of the boundary conditions, the numerical
resolution, the hydrodynamics scheme, the tensor viscosity, the choice of ray
directions for the radiative transport, and the time-step control (Freytag et
al., 2012; Freytag, 2013). One particular point concerns the closed-bottom
boundary: as the radiative flux in the convection zone of cool stars with
efficient convection is negligible, the energy in the bottom layers is not
injected by radiation but is added as an energy source term within the central
“internal” boundary. In these layers, the presence of strong velocity fields
(waves or downdrafts) needs an adapted treatment of boundary conditions to
avoid the generation of spurious acoustic waves when downdrafts bounce against
a hard lower boundary. In addition to this, there are also options to set the
boundary conditions handling the energy-source terms with the possibility to
include a drag force in the core (Freytag, 2017). The (magneto-) hydrodynamics
and radiation-transport solvers integrate right through this zone.

Figure 1: Example of a starting RHD simulation for a RSG star (Chiavassa et
al., 2009). The panels show the evolution of the bolometric intensity from a
sphere in hydrostatic equilibrium (leftmost panel) to the relaxed simulation
(rightmost panel). The artifacts caused by the mismatch between the spherical
object and the Cartesian grid become less evident with time.

The initial model is a sphere in hydrostatic equilibrium with a weak velocity
field inherited from a previous model with different stellar parameters. The
initial temperature stratification follows a gray atmosphere with the
optically thick layers chosen to be adiabatic. Fig. 1 shows the time
evolution. The first snapshot in the leftmost panel does not display any
convection-related structures but with some regular patterns due to the coarse
numerical grid and poor temperature sampling at the sharp edge of the
atmosphere. Later, a central spot appears and then vanishes completely when
convection becomes strong. After several stellar years, the intensity contrast
has grown and the surface pattern becomes completely irregular. All memory
from the initial symmetry is lost as well as the influence of the Cartesian
grid. Once the simulation has reached a relaxed state, stellar parameter
fluctuations stabilize and average quantities can be computed.

Table 1 displays the temporal and spherically averaged stellar parameters and
the numerical box details of two RHD simulations representing AGB and RSG
stars, while Fig. 2 shows the spherical averaged quantities still vary as a
function of time as explained in Chiavassa et al. (2009) and Chiavassa et al.
(2011a).

---
Figure 2: Spherically averaged, global properties as a function of time for an
AGB (Left-hand column) and an RSG simulation (Right-hand column). See Table 1
for their temporal averages. Table 1: Parameters of the RHD simulations. The
first two columns show the simulation name and stellar type. The following
five columns display the stellar parameters such as the mass ($M_{\star}$)
used for the external potential, the average emitted luminosity ($L_{\star}$),
the Rosseland radius ($R_{\star}$), the effective temperature
($T_{\mathrm{eff}}$) and surface gravity ($\log g$) at the Rosseland radius.
Then the pulsation period $P_{\mathrm{puls}}$ and the full-width half-maximum
(FWHM) of the distribution of pulsation frequencies $\sigma_{\mathrm{puls}}$.
All these quantities are averaged over different epochs (9th column,
$t_{\mathrm{avg}}$). Errors are one-standard-deviation fluctuations with
respect to the time average (see Chiavassa et al., 2011a, 2009). Solar
metallicity is assumed. The simulation of the AGB comes from the grid of
Freytag et al. (2017) and it is used also in Chiavassa et al. (2021). The
simulation of the RSG has been computed including a rotational period of 30
years adding a centrifugal potential to the gravitational potential and the
Coriolis force which is applied at each hydrodynamics step (Freytag et al.,
2017).

Simulation | $M_{\star}$ | $L_{\star}$ | $R_{\star}$ | $T_{\mathrm{eff}}$ | $\log g$ | $P_{\mathrm{puls}}$ | $\sigma_{\mathrm{puls}}$ | $t_{\mathrm{avg}}$ | Grid | xbox
---|---|---|---|---|---|---|---|---|---|---
| $[M_{\odot}]$ | $[L_{\odot}]$ | $[R_{\odot}]$ | [K] | [cgs] | [yr] | [yr] | [yr] | [grid points] | $[R_{\odot}]$
st29gm06n001 | 1 | 6956.3$\pm$547.4 | 350.4$\pm$11.3 | 2814.7$\pm$68.9 |  $-$0.65$\pm$0.03 | 1.150 | 0.314 | 25.35 | 2813 | 1381
AGB |  |  |  |  |  |  |  |  |  |
st35gm04n43 | 5 | 41033.4$\pm$1850.1 | 616.2$\pm$8.4 | 3308.2$\pm$33.2 |  $-$0.45$\pm$0.01 | 1.456 | 0.209 | 30.21 | 3153 | 1626
RSG |  |  |  |  |  |  |  |  |  |

####  3.1.2 Global quantities

---
Figure 3: Logarithm of temperature (_left column_), density (_central
column_), and opacity (_right column_) of central slice of a snapshot for the
RSG (left column) and AGB (right column) simulations of Table 1.

RHD simulations of evolved stars show a very heterogeneous surface caused by
the dynamical granulation that evolves on timescales of weeks to years
(Chiavassa et al., 2011b; Freytag et al., 2017; Ahmad et al., 2023). The
radiation is of primary importance for many aspects of convection and envelope
structure. Radiation cools down the surface to provide a somewhat unsharp
outer boundary for the convective heat transport and it also contributes
significantly to the energy transport in the interior, where convection never
carries close to 100% of the luminosity (Chiavassa et al., 2011a). Moreover,
it must be noted that the entropy jump, happening in the optically thick
region from the atmosphere to the layers below, is fairly large.
Below the visible layers (i.e., optical depth $\tau_{\mathrm{Rosselend}}>1$),
the opacity has its peak at $T\sim 13,000$ K (see the temperature and opacity
surface rendering in Fig. 3) causing a very steep temperature jump which is
prominent on top of upflow regions. At the same time a density inversion
appears (see density in Fig. 3), which is a sufficient condition for
convective instability. The entropy drop occurs in a very thin layer, while
the smearing from the averaging procedure over nonstationary up- and downflows
leads to the large apparent extent of the mean superadiabatic layer (Chiavassa
et al., 2011a).

Figure 4: Logarithm of the density (increasing from blue to red) in a slice
through the center of the RSG simulation of Table 1 overplotted to the
isosurface of the Rosseland photosphere ($\tau_{\rm Ross}=1$), where the
continuum flux is formed, is shown in amaranth color.

The emerging visible luminosity is related to layers where waves and shocks
dominate together with the variation in opacity through the atmosphere. The
global convective flow in the simulations is characterised by rising material
that originates in the NBconvective zone below the Rosseland radius and
develops as an atmospheric shock when it reaches higher values in radius
(Freytag et al., 2017; Liljegren et al., 2018). In deeper layers, non-
stationary convection (e.g., merging down-drafts or other localized events)
produces sound wave in the stellar interior and only a few large-scale cells
dominate the convective flow. The sound waves travel through the stellar
interior ($\tau_{\mathrm{Rosselend}}>1$) to the outer optically thin layers
where they hit the surface, slow down and get compressed due to the drop in
temperature and sound speed. The amplitude of the sound wave rises due to the
decrease in density (Fig. 4) and turns into a shock (Kravchenko et al., 2019).
Here, the turbulent pressure can be as high as more than 20 times the gas
pressure (Fig. 5, top row). Just above the stellar radius, the turbulent
pressure contributes significantly to the average pressure stratification and
to the radial (supersonic) velocities.

---
Figure 5: _Top row:_ ratio between turbulent pressure and gas pressure for
different snapshots of the simulations. _Bottom row:_ spherical shells
averages over time of the radial Mach numbers. The left panels show the RSG
and the right one the AGB simulation. The red vertical lines in all the panels
is the approximative position of the radius from Table 1.

The newly formed shocks propagate from the stellar surface and out of the
simulation domain with significant mach numbers of 3–4 on average (Fig. 5,
bottom row). These shock waves have lasting effects on the outer layers making
them extremely heterogeneous and having Mach numbers peaking at 8–10. (Fig.
6). In these layers both the density and the temperature show irregular
structures with convection cells in the interior and a network of shocks above
the photosphere (Fig. 4 and Fig. 3). Local fluctuations in high Mach numbers
as well as small-scale heights due to shocks pose high demands on the
stability for the hydrodynamic solver. A side effect of the steep and
significant temperature jump is the increase in pressure scale height from
small atmospheric values to values that are a considerable fraction of the
radius in layers just below the atmosphere. The simulation in Fig. 6 shows
that the outer boundaries are either hit at some angle by an outgoing shock
wave or let the material fall back (mostly with supersonic velocities). The
shocks pass through the boundaries with a simple and stable prescription in
the code based on filling typically two layers of ghost cells where the
velocity components and the internal energy per unit mass are kept constant
(Freytag et al., 2012).

Figure 6: Radial Mach numbers computed for the RHD simulation of the AGB in
Table 1 (blue to red colors) rendered transparently over the isosurface of the
temperature equal to 2800 K (shaded, gray surface). This temperature
corresponds approximatively to the effective temperature of the simulation
(Table 1).

###  3.2 The convective envelope as seen by Athena++

In this section, we focus on the Athena++ simulations which have been used to
calibrate mixing length parameters and predict signatures of the 3D convective
outer layers in the supernova shock breakout (Section 4.2). We describe the
stellar structure as seen in the simulations of two RSG envelopes (Table 2) in
higher- and lower- luminosity regimes.

####  3.2.1 Conducting a simulation in Athena++

Name | $R_{\mathrm{IB}}/R_{\odot}$ | $R_{\mathrm{out}}/R_{\odot}$ | heat source | resolution ($r\times\theta\times\phi$) | duration | $m_{\mathrm{IB}}/M_{\odot}$ | $M_{\mathrm{final}}/M_{\odot}$
---|---|---|---|---|---|---|---
RSG1L4.5 | 400 | 22400 | “hot plate” | $384\times 128\times 256$ | 5865 d | 12.8 | $16.4$
RSG2L4.9 | 300 | 6700 | fixed $L$ | $256\times 128\times 256$ | 5766 d | $10.79$ | $12.9$

Table 2: Simulation properties of the Athena++ models from Goldberg et al.
(2022a), including inner boundary radius ($R_{\mathrm{IB}}$), outer boundary
($R_{\mathrm{IB}}$), heat source, resolution, run duration, core mass
$m_{\mathrm{IB}}$, and total mass at the simulation end
($M_{\mathrm{final}}$). All models have $\theta=\pi/4-3\pi/4$ and
$\phi=0-\pi$, with $\delta r/r\approx 0.01$. The naming scheme indicates the
log of the time-averaged luminosity $\log(L/L_{\odot})$.

The Athena++ Red Supergiant setup described in Goldberg et al. (2022a) is
based on the semi-global (large solid-angle) “Star3D” setup presented in Jiang
et al. (2018), with the difference that the gravitational potential accurately
includes the (non-negligible) time-dependent, spherically symmetric, envelope
mass, $g(r)=-\nabla\Phi(r)=G(m_{\rm IB}+m_{\rm env}(r))/r^{2}$ where $m_{\rm
env}(r)$ is the mass interior to the shell at radius $r$. These simulations
cover 70.6% of the face-on hemisphere (i.e. solid angle $\Omega=1.41\pi$)
using spherical polar coordinates with 256 uniform bins in azimuth $\phi$ from
$0-\pi$ and 128 uniform bins in polar angle $\theta$ from $\pi/4-3\pi/4$, with
periodic boundary conditions in $\theta$ and $\phi$ enforced by ghost zones
outside the simulation domain (see Stone et al., 2020) to conserve mass and
energy flux and to avoid the “splashback” phenomenon associated with
convective motions through reflective boundaries, and logarithmic bins in
radius with $\delta r/r\approx 0.01$.

Figure 7: History of the averaged radial profiles for the RSG1L4.5 (red, left)
and RSG2L4.9 (orange, right) Athena++ models. Top to bottom show log(density),
magnitude of the 3D velocity, log(temperature), and log(opacity). All
logarithms are base 10 and units inside the logarithms are cgs, and velocites
are reported in km/s. Vertical dotted lines indicate when the envelope appears
to have reached a convective steady state. Horizontal dashed lines approximate
the region where some fraction of the stellar area has $\tau>\tau_{\rm crit}$.
Figure from Goldberg et al. (2022a).

The inner boundary (IB) has fixed radius at $400R_{\odot}$ and $300R_{\odot}$
for RSG1L4.5 and RSG2L4.9 (named by the log of the time-average luminosity
leaving the simulation domain) respectively with a zero-velocity, fixed
temperature, fixed density inner boundary condition, and analysis is
restricted to discussing the convective region which exhibits steady-state
behavior, from $\approx$ 450$R_{\odot}$ outwards. The implicit radiation
scheme implemented in Athena++ (Jiang et al., 2014; Jiang, 2021) is used,
which directly solves the RHD equations for specific intensity over discrete
ordinates, including radiation in both the energy and momentum equations as a
source term. The gas equation of state is taken to be a gamma-law where
$\gamma=5/3$, as the radiation pressure is included explicitly. For the two
models discussed, initial conditions are constructed in radiative equilibrium
from these values at the inner boundaries which were motivated by 1D models
constructed using the stellar evolution code MESA (Paxton et al., 2011, 2013,
2015, 2018, 2019; Jermyn et al., 2023). The luminosity is then either supplied
by a “hot plate” method (RSG1L4.5), where the temperature is increased at the
inner boundary to drive convection, or a fixed high L coming in at the inner
boundary (RSG2L4.9).

The simulations are run until they achieve a convective steady state as
discussed in more detail in Goldberg et al. 2022a. This time is shown by the
vertical dotted line in Fig. 7, which presents the time-evolution of the
radially-averaged density, velocity, temperature, and opacity for the two
Athena++ models. Also indicated by the horizontal dashed line is the location
where radiation is expected to carry a significant fraction of the stellar
flux, $\tau\approx\tau_{\rm crit}$ where $\tau_{\rm crit}\approx P_{\rm
rad}/P_{\rm gas}c/|v_{r}|~{}$ is a few hundred in the RSG regime (e.g.
Goldberg et al. 2022a; Jermyn et al. 2022; Schultz et al. 2023).

####  3.2.2 Envelope Structure and Convective Properties in Athena++

Like the CO5BOLD models discussed above, the Athena++ RHD simulations show
very heterogeneous outer layers with convective overturn timescales on the
order of a year.

Figure 8:  Equatorial ($\theta=\pi/2,z=0$) slice of a characteristic snapshot
of the RSG2L4.9 model, showing (upper half of the figure:) radial velocities
and (lower half of the figure:) density fluctuations relative to the shellular
volume-averaged density at each radius. The simulation domain is from $\phi=0$
to $\pi$; thus the image is reflected about $y=0$ as indicated by the axis
labels. Figure from Goldberg et al. (2022a). Figure 9: Left: Surface rendering
of the radial velocity fluctuations at $r=600R_{\odot}$ in RSG1L4.5 at day
4707. Right: Snapshot of radial velocity slices for the same model. Panels
show radial slices at $r=500R_{\odot}$ (top) to $800R_{\odot}$ (bottom) in
$100R_{\odot}$ intervals, and axes show the extent in azimuth $\phi$ and co-
polar angle cos($\theta$). The volume-weighted mean sound speed is 33 km/s at
$r=500R_{\odot}$, 26 km/s at $600R_{\odot}$, 19 km/s at $700R_{\odot}$, and 8
km/s at $800R_{\odot}$. Figure adapted from Goldberg et al. (2022a).

The radial and convective velocities in the Athena++ models remain coherent on
large scales throughout the convective envelope (see Fig. 9), in broad
agreement with the flow morphology of the CO5BOLD simulations (e.g., recent
work by Ma et al. 2023). The density fluctuations are $\approx 10\%$ in the
middle of the convective envelope, reaching unity and greater at fixed radius
in the outer envelope, as seen in Fig. 8 for RSG2L4.9. Although the net
angular momentum in the envelope is nearly zero, at any given time the large-
scale tangential velocity fluctuations (tens of km/s, coherent across tens to
hundreds of solar radii) result in finite specific angular momentum
($j_{\mathrm{rand}}\sim 10^{18}-2\times 10^{19}$cm2/s) at each radius,
consistent with the recent 3D hydrodynamic simulations of interior RSG
convection Antoni and Quataert (2022).

Fig. 10 shows profiles of the entropy, radiation flux, and kinetic energy to
thermal energy ratio in the Athena++ RSG models of Goldberg et al. (2022a). As
radiation is included in both the energy and momentum equations in the
Athena++ setup, the simulations achieve a flat entropy profile in the
interior, a signature of efficient convection (as predicted by 1D models and
MLT) which decreases and steepens as the convection becomes increasingly
super-adiabatic. Radiation only carries $10-20\%$ of the flux in the interior,
with convection carrying the rest. In the near-atmospheric layers, turbulent
mach numbers in the vigorous, inefficient convection cross unity and the
kinetic energy density matches or exceeds the thermal energy content, with
turbulent pressure contributing significantly to hydrostatic balance.

Figure 10:  Specific entropy (top panels), radiative luminosity (middle
panels), and ratio of the turbulent kinetic energy to the internal energy
(bottom panels) in the RSG1L4.5 model at day 4707 (red, left) and RSG2L4.9 at
day 4927 (orange, right). Mass-weighted averages are shown in black, with 80%
of the mass lying within the shaded regions. The 1D photosphere, where
$\langle{L(r)}\rangle=4\pi
r^{2}\sigma_{\mathrm{SB}}\langle{T_{r}(r)}\rangle^{4}$, is given by the
vertical dashed line. Figure from Goldberg et al. (2022a).

The fluctuations in the stellar properties, which increase radially outwards
in the envelope, can be seen in Fig. 11. The ratio of radiation to gas
pressure is also significant in the interior, with $\langle P_{\rm rad}/P_{\rm
gas}\rangle\approx 0.15-0.5$ at the $r=450R_{\odot}$ coordinate, decreasing to
the level of a few percent approaching the photsphere. Rather than a smooth
transition from the H opacity peak to scattering and neutral opacities, the
temperature and opacity show bimodal behavior (not present in the density
profiles) in the $\approx 100-200R_{\odot}$ region beneath the photosphere for
both models. Thus the sharp opacity drop seen in 1D models which use only the
volume-averaged $\rho$ and $T$ profiles appears less steep in the 3D
simulations, though it is equally sharp along any given line-of-sight, as seen
also in the CO5BOLD simulations.

Figure 11: Top to bottom: Density, temperature, $P_{\rm rad}/P_{\rm gas}$
ratio, and opacity for the RSG1L4.5 model at day 4707 (red, left panels) and
the RSG2L4.9 model at day 4927 (orange, right panels). Color saturation
indicates the volume-weighted probability of finding a fluid element at a
given ($\rho,T,$ pressure scale height $H/r,P_{\rm rad}/P_{\rm gas},\kappa$)
at each radial coordinate. Solid black lines give the volume-weighted averages
of each (non-log) quantity ($\langle\rho\rangle,\langle T\rangle,\langle
P_{\rm rad}/P_{\rm gas}/\rangle,\langle\kappa\rangle$). The vertical black
dashed line is the photosphere radius $R_{\rm phot}$. The ratio $P_{\rm
rad}/P_{\rm gas}$ is only shown for $r\leq R_{\rm phot}$. The $\kappa$ panels
(bottom) show both the volume-averaged opacity $\log(\langle\kappa\rangle)$
reported by Athena (solid lines), as well as the OPAL opacity (dotted lines)
recovered by inputting the volume-averaged $T$ and $\rho$ profiles. For
reference, the cyan line in the upper left panel shows the density profile of
the MESA model which motivates the inner boundary and simulation initial
conditions. Figure from (Goldberg et al., 2022a).

Notably, the opacity $\kappa_{\mathrm{Edd}}$ above which $L$ locally exceeds
the Eddington Luminosity $L_{\rm Edd}$, $\kappa_{\mathrm{Edd}}=4\pi GcM/L$, is
$\approx 1-10$ cm2/g, which means that a significant amount of the material at
those transitionary radii with $\kappa\sim$ tens of cm2/g
$>\kappa_{\mathrm{Edd}}$. At optical depths with radiation-dominated energy
transport where some material exists at $\tau_{\rm crit}>\tau>1$, radiation
forces may significantly impact fluid motion with high
$\kappa\gtrsim\kappa_{\rm Edd}$, and MLT-like convection is not expected. In
fact in those layers of the Athena++ simulations, the typical correlations
between the radial velocity and the density, opacity, and entropy flip; such
that the cold (low-entropy), opaque, dense regions are correlated with outward
velocities, rather than simply sinking like they would in MLT-like convection.
Similar inverted-correlation behavior is also seen in other simulations of
luminous stars (e.g. in OB-star envelopes; Schultz et al. (2022)), but not in
simulations of solar-like convection (e.g. Stein and Nordlund 1998), likely
owing to RHD effects where $\tau_{\rm crit}\gg 1$ and $L\gtrsim L_{\rm Edd}$.

##  4 Results and applications

###  4.1 CO5BOLD applications: Convection cycles in evolved stars

Freytag et al. (2017) and Ahmad et al. (2023) reported that convection and
pulsations are both emergent properties of the RHD simulations. Moreover,
Ahmad et al. (2023) showed the non-linear nature of the interplay between
them: in all of their 3D RHD simulations, the radial pulsations can be
observed and the pulsation periods in the fundamental mode extracted, with
values in agreement with the observed periods of long period variable (LPV)
stars. The authors also highlighted the importance of the characteristics of
deep convective regions and large convection cells, for understanding the
interplay with pulsations.
In this Section, we concentrate on the convection cycles based on the
observational approaches related to visual magnitude and radial velocities
measurements. Measurements of radial velocities are fundamental to determining
stellar space velocities to investigate the kinematic structure of stellar
populations in the Galaxy or to monitor for radial velocity variations, either
of which could point to the presence of an unseen companions. In addition to
this, radial velocities can also be used to recover the projected velocity
field at different optical depths in the stellar atmosphere using the
tomographic method described in the Section 4.1.2.

In this context, convection plays an important role in the formation of
observed spectral lines and deeply influences the spectral line formation in
late-type stars (Asplund et al., 2000). Each line, forming across a range of
depths in the atmosphere, has a unique fingerprint in the spectrum that
depends on line strength, depth, shift, width, and asymmetry across the
granulation pattern depending on their height of formation and sensitivity to
the atmospheric conditions of main sequence and up to red giant branch stars
(Allende Prieto et al., 2013; Chiavassa et al., 2018a). The RHD simulations
provide a self-consistent ab-initio description of the non-thermal velocity
field generated by convection, shock waves, and overshoot that manifests
itself in spectral line shifts and changes in the equivalent width (Chiavassa
et al., 2011a; Chiavassa and Freytag, 2015) as revealed in the observations
(e.g. Gray 2008).

Finally, the velocity field in the evolved stars, reconstructed using the
tomographic method (Section 4.1.2), can be related to the photometric
variability to reveal phase shifts resulting in a hysteresis loop (Section
4.1.3. According to Gray (2008), the hysteresis loop illustrates the
convective turn-over of material in the stellar atmosphere: first, the rising
hot matter reaches upper atmospheric layers, then the temperature drops as the
matter falls back into the star.
In the following part of this Section, we present a detailed discussion of the
methodology and possible origin of the hysteresis loop in evolved cool stars,
taking into account the mechanisms related to the non-stationary effects of
convection.

####  4.1.1 Pulsation properties

The irregular structures of convection cells in the interior are visible in
the temperature and density structures (Fig. 3) as well as in the radial
velocity maps (Fig. 12 and 13) together with a network of shocks in the
atmosphere of both the RSG or the AGB simulations. The high and heterogeneous
velocities (up to $\sim$35 km/s in RSG simulation) and pressure fluctuations
determine the acceleration of shock waves and, eventually, affect the mass
loss. For the RSG, the large shocks just above $\tau_{\mathrm{Rosselend}}$ =1
can be $\sim 200-250\,R_{\odot}$ wide (each grid point is about
$5\,R_{\odot}$, Table 1) for both the RSG and AGB with a local temperature of
$\sim$2500 K and $\log(\rho/[{\rm g\,cm}^{-3}])\sim-13$. This width is
comparable to the actual radius of the star.
The rising material (blue color) observed in the figure originates in the deep
convective zone and an atmospheric shock develops higher up
($\tau_{\mathrm{Rosselend}}$ lower than 1) along the following steps. Freytag
et al. (2017) and Liljegren et al. (2018) illustrated that in AGB simulations
(including the one of Table 1) non-stationary convection (e.g., merging down-
drafts or other localized events) manifests into giant convective cells, which
change topology on very long timescales, and short-lived small surface
granules. Moreover, the authors managed to measure the cycles of outward
moving shocks and material falling back toward the star and quantify a period
for the pulsation modes.

The radial motions in the atmosphere can then be retrieved using averages over
spherical shells of the radial velocities for each snapshot of the simulation
(Fig. 14, top panels) for the RSG simulation of Table 1. As in Freytag et al.
(2017) and Ahmad et al. (2023), the behavior of the inner part of the model
differs from that of the outer layers: below $\sim 600\,R_{\odot}$ (the
nominal radius is $616\,R_{\odot}$) the velocity field is rather regular and
coherent over all layers, close to a standing wave. The different slopes
visible in the outer layers (above $\sim 600\,R_{\odot}$) are clearly
indicating the presence of propagating shock waves but in a less regular and
smooth way than in the AGB case (as a matter of comparison, see Fig. 5 of
Freytag et al. 2017). The periodic behavior of this movement can be extracted
using a Fourier analysis and plotted (Fig. 14, bottom panel). A clear
prominent frequency is noticeable at $R\sim 600\,R_{\odot}$, while far outside
in the shock-dominated atmosphere the possible other modes are significantly
smaller in amplitude than in the interior of the star. A Gaussian distribution
was fitted in the frequency domain to characterize the largest signal. The
resulting period of pulsation is $P_{\mathrm{puls}}=1.456\pm
0.209\,\mathrm{yr}$.
These radial pulsations are not transient phenomena introduced by the initial
conditions but the outcome of a mode-excitation mechanism caused by the
convective motions in the stellar interior, as confirmed for the first time
for AGB simulations by Freytag et al. (2017).

---
Figure 12: _Top panel:_ slice across the numerical box showing the radial
velocity for the RSG simulation of Table 1. The blue indicates outward and red
inward flow. Moreover, the temperature volume rendering is also displayed
(yellow to black colors). The temperature values arbitrarily stop at the
effective temperature of the simulations to show the approximate position of
the $\tau_{\mathrm{Rosselend}}=1$. In the simulations the temperature can go
down to $\sim$1000 K. _Bottom panel:_ Enlargement of the radial velocity slice
from above.
---
Figure 13: Same as in Fig. 12 for the AGB simulation of Table 1.
---
Figure 14: _Top panels:_ spherically averaged radial velocities for the whole
run time and radial extent of the RSG simulation of Table 1. The different
colors show the average vertical velocity at that time and radial distance.
The rightmost panel displays the area in the dashed rectangle of left panel
with mass-shell movements plotted as iso-mass contour lines. The velocity
range and color is the same as in Fig. 12. _Bottom panel:_ Normalised power
spectra derived from the velocity fields from above, mapped over frequency and
radial distance.

####  4.1.2 The tomographic method

The study of the velocity field in evolved stellar atmospheres poses a special
challenge, as the spectra of these stars are extremely crowded, particularly
in the optical domain. Moreover the combined convective and pulsational
velocity field have a large impact on the formation of spectral lines
producing asymmetries and Doppler shifts (Fig. 15). To disentangle various
contributors to the complex line shape, high-resolution spectroscopic
measurements coupled with the cross-correlation technique provide a powerful
tool to overcome these difficulties. The information relating to the line
doubling (i.e. two distinct velocity components, contributing to the line
shape) is in fact distributed among a large number of spectral lines, and can
be summed up into an average profile, or more precisely into a cross-
correlation function (CCF). If the end correlation of the stellar spectrum
with a mask involves many lines, it is possible to extract the relevant
information from very crowded and/or low signal-to-noise spectra. Moreover,
radial velocities must be traced on both spatial and temporal scales to allow
detailed reconstruction of atmospheric motions in evolved stars. A recently
established tomographic method is an ideal technique for this purpose. It must
be stated right away that the word ”tomography” is used here with a meaning
close to its etymological roots (”slice, section”).

---
Figure 15: Example of a synthetic spectral line of the Ti I at 6261.11 Å for a
snapshot of a 3D RSG simulation (Chiavassa and Freytag, 2015). The vertical
dashed line shows the spanned velocities of the line bisector. The different
arrows and colors display the approximate position of different velocity
components, resulting from the heterogeneous velocity field (Fig. 12).

The tomographic method (Alvarez et al., 2001) aims to recover the temporarily-
resolved line-of-sight velocity distribution as a function of depth in the
stellar atmosphere. The technique is based on sorting spectral lines according
to their formation depths. The atmosphere is split into different slices
corresponding to specific optical depth ranges. For each slice, a spectral
mask is then constructed, which contains the wavelengths of lines forming in
the corresponding range of optical depths. The cross-correlation functions
(CCFs) of the masks with stellar spectra reflects the average shape and radial
velocity (RV) of lines forming in a given optical depth range.

The concept of the tomographic method originates from the recognition by
Schwarzschild (1952) that the doubling of the H-$\alpha$ line observed in
Cepheid variables of the W Vir type follows a specific time sequence that may
be explained by a shock front passing through the line-formation zone. Line
doubling has since been observed in many other lines, both in the infrared and
in the optical. For long-period variables (LPVs) of Mira type, Merrill (1955)
was the first to suggest that the observed spectral changes may be explained
by the presence of a shock wave moving outwards. In pulsating atmospheres of
LPVs or Cepheids, what has since been known as the Schwarzschild scenario
(Schwarzschild, 1952, see Fig. 16) states that, when the shock front is
located far below the line-forming region, the stellar atmosphere contains
only falling material so that the spectral lines are red-shifted. Then, as the
shock wave lift up, a blue component appears and strengthens. This outward
motion of the shock front has been revealed by the tomographic method applied
to Mira variables (Fig. 17). One example is reported in Alvarez et al. (2000)
for a particular AGB stars, namely Mira, for which the tomographic method
clearly unveils the Schwarzschild scenario as a line-doubling (ie., the shock
signature, Fig. 16) appearing in diverse layers of time-consecutive CCF
probing, from below to above the stellar surface. In other words, the
tomographic method revealed the upward motion of the shock fronts combining
temporal and spatial variations of the CCFs (Jorissen et al., 2016).

The tomographic method of Alvarez et al. (2001) was applied to short time-
series spectra of RSG stars by Josselin and Plez (2007). This study detected
line asymmetries in resulting cross-correlation profiles and the absence of
the Schwarzschild scenario (Fig. 18). Thus, RSGs are interesting targets for a
more detailed tomographic study.

Figure 16: Illustration of the Schwarzschild scenario: temporal sequence
following the intensity of the red and blue components of absorption lines,
when a shock wave propagates through the photosphere. Figure from Alvarez et
al. (2000). Figure 17: Sequence of cross-correlation profiles of the Mira star
RT Cyg obtained with four tomographic masks, around maximum light. Columns
correspond to different masks, the leftmost mask probes the innermost
atmospheric layer. Rows correspond to variability phases varying from 0. to 1.
The cross-correlation profiles follow the Schwarzschild scenario: a single red
component is progressively transforming into a single blue component.
Reproduced from Jorissen et al. (2016).  Figure 18: Sequence of cross-
correlation profiles of the RSG stars SW Cep (left) and $\mu$ Cep (right)
obtained for eight tomographic masks. Rows correspond to different masks, mask
C1 probes the innermost atmospheric layer. Columns correspond to different
observing epochs. In contrast to Mira stars, the CCFs of RSGs do not follow
the Schwarzschild scenario. Figure from Josselin and Plez (2007).

A detailed description of the tomographic technique and its recent
developments can be found in Kravchenko et al. (2018). The method relies on
the design of spectral templates (so-called masks) containing lines forming at
given, pre-specified ranges of optical depths. The original implementation of
the method (Alvarez et al., 2001) relied on the Eddington-Barbier
approximation, stating that in a grey atmosphere, the flux emerges from the
layer located at optical depth $\tau_{\lambda}=2/3$. Thus a line at wavelength
$\lambda$ was supposed to form at optical depth $\tau_{\lambda}=2/3$. An
appropriate model atmosphere is then used to convert this monochromatic
optical depth $\tau_{\lambda}=2/3$ into a reference optical depth (at 5000
$\AA$ for instance), which is then a proxy for geometrical depth. As shown by
Magain (1986), the Eddington-Barbier approximation gives however a valid depth
of line formation only for sufﬁciently strong lines. In the recent
implementation of the tomographic method (Kravchenko et al., 2018), the
contribution function to the line depression is used in order to correctly
assess the depth of formation of spectral lines. Albrow and Cottrell (1996)
derived the contribution function111The amount of light contributed to the
surface flux as a function of depth in the atmosphere (CF) for a single ray
(SR) of a 1D plane-parallel atmosphere:

${\rm CF_{SR}}(\log\tau_{0},{\lambda})=(\ln
10)\,\mu^{-1}\frac{\tau_{0}}{\kappa_{c,0}}\kappa_{l_{\lambda}}\,(I_{c_{\lambda}}-S_{l_{\lambda}})\,e^{-\tau_{\lambda}/\mu}$
(1)

and for the whole stellar disk (DI):

${\rm CF_{DI}}(\log\tau_{0},{\lambda})={2\pi}\,(\ln
10)\frac{\tau_{0}}{\kappa_{c,0}}\int_{0}^{1}\kappa_{l_{\lambda}}\,(I_{c_{\lambda}}-S_{l_{\lambda}})\,e^{-\tau_{\lambda}/{\mu}}\,d\mu$
(2)

In the above equations, $\tau_{0}=\int\kappa_{c,0}\,\rho\,{\rm d}x$ is the
continuum optical depth at a reference wavelength $\lambda_{0}=5000$ Å, $\rho$
is the density, $\kappa_{c,0}$ is the continuum absorption coefficient at
$\lambda_{0}$, $I_{c}$ is the continuous intensity, $I_{l}$ is the line
intensity, $\kappa_{l}$ is the line absorption coefficient, $S_{l}$ is the
line source function, $\tau_{\lambda}=\int\kappa_{\lambda}\,\rho\,{\rm d}x$ is
the optical depth along the ray, $\kappa_{\lambda}$ is the total absorption
coefficient at wavelength $\lambda$, $\mu=\cos{\theta}$ is the cosine of the
angle $\theta$ between the line of sight and the radial direction ($\mu=1$ for
the direction towards the observer). The $\tau_{0}$ scale can be used as a
proxy to the geometrical depth $x$.

Kravchenko et al. (2018) implemented the computation of the CF in the 1D
radiative transfer code TURBOSPECTRUM (Plez, 2012), which performs spectrum
synthesis from static one-dimensional (1D) MARCS model atmospheres (Gustafsson
et al., 2008). Fig. 19 illustrates an example of the CF emerging from a 1D
hydrostatic model of RSG star. The set of $\tau_{0}$ corresponding to the
maximum CF for each wavelength defines a crest line called ”depth function”.
The depth function resembles the shape of a spectrum but instead of intensity
(or flux) it indicates the central optical depths $\tau_{0}$ at which the
emergent intensity at wavelengths $\lambda$ forms (left panel of Fig. 20).
This allows to split the atmosphere into different optical-depth layers
(colored bands in Fig. 20) and construct tomographic masks across the stellar
atmosphere. Each mask contains the central wavelengths of spectral lines whose
depth function minimum falls within a range of optical depths defined by a
given layer. The cross-correlation of the masks with observed or synthetic
stellar spectra provides the distribution of the line-of-sight velocity field
inside a given atmospheric slice (Fig. 21).

---
Figure 19: An example of the Contribution Function emerging from a 1D
hydrostatic model of a RSG star (figure from Kravchenko et al., 2018). The
crest line shown in red defines the depth function illustrated in Fig. 20.
---
Figure 20: Sketch of the mask construction process. _Left panel:_ The depth
function (top) and the normalised synthetic flux (bottom) are shown as a
function of wavelength for a 1D RSG model. The depth function corresponds to
the crest line of the CF shown in Fig. 19. The comparison of the depth
function with the synthetic spectrum shows that the cores of spectral lines
form in outer layers while the wings form in deeper layers. The different
colors sketch the approximate position of the masks used to divide the
atmosphere. _Right panel:_ The different atmospheric layers spanned by the
masks are sketched over a slice of a RSG RHD simulation representing the line-
of-sight velocity ($V_{z}$). Note that the masks are not assumed to be
concentric spherical shells. Figure 21: Example of a temporal sequence
(columns) of cross correlation functions (CCF, obtained by applying a set of
five tomographic masks to a series of observed stellar spectra. The masks (C1
to C5) have the same colors as in Fig. 20.

####  4.1.3 Recovering velocity fields across the atmosphere in RSG stars

The validation of the tomographic method on 3D RHD simulations needs the post-
processing calculations of high spectral resolution synthetic spectra from the
3D thermodynamical structures. For this purpose, a possible option is the use
of the 3D pure-LTE radiative transfer code Optim3D (Chiavassa et al., 2009).
This code computes synthetic spectra and intensity maps for any snapshots of
RHD simulations and uses pre-tabulated extinction coefficients generated with
the MARCS code (Gustafsson et al., 2008) and by adopting the solar composition
of, e.g., Asplund et al. (2009). Then, the extinction coefficients are
interpolated at each Doppler-shifted wavelength in each cell along the ray.
The numerical integration of the equation of radiative transfer gives the
intensity emerging towards the observer at that wavelength and position. This
calculation is performed for every $x$-/$y$-point of a face of the simulation
cube, and for all the required wavelengths. Using these spectra, the
tomographic method has been tested by recovering the velocity field
distribution across the atmosphere of a 3D RHD simulation of an RSG star in
Kravchenko et al. (2018). The authors derived the 3D CF for the ”star-in-the-
box” simulation, which was implemented in the radiative transfer code Optim3D:

$CF_{\rm 3D}(\tau_{0})={\rm
ln}10\sum_{i,j=0}^{N}\,\tau_{0}\,\frac{\kappa_{l,i,j}}{\kappa_{0,i,j}}\,\bigg{[}I_{c,z}\Big{(}x_{i},y_{j},z(\tau_{0})\Big{)}-\\\
S_{l,z}\Big{(}x_{i},y_{j},z(\tau_{0})\Big{)}\bigg{]}e^{-\tau}\frac{\Delta
x_{i}\Delta y_{j}}{R^{2}}$ (3)

where $N$ is the number of grid points in the numerical box along the line-of-
sight $z$-axis, $i$ and $j$ are ray numbers along the x and y axes. For the
equidistant grid $\Delta x_{i}=\Delta y_{j}$. All quantities in the equation
depend on wavelength $\lambda$. Fig. 22 illustrates an example of the 3D CF
computed with Eq. 3 for the same wavelength range as the 1D CF from Fig. 19.
Similar to the 1D CF, the 3D CF is characterized by a single maximum in the
$\log\tau_{0}$ direction.

---
Figure 22: An example of the 3D CF for a 3D RHD simulation of a RSG star from
Table 1 for the same wavelength range as Fig. 19 (figure from Kravchenko et
al., 2018).
---
Figure 23: Left panel: Distribution of the line-formation depths in a 3D RSG
snapshot (blue color) for a set of eight tomographic masks. The red bands
represent the optical depth ranges used for selecting lines in a given mask
from the 1D contribution function. Right panel: Reconstruction of the velocity
field in the 3D snapshot as a function of the optical depth in the atmosphere.
Adopted from Kravchenko et al. (2018).

To test the tomographic method, a set of eight tomographic masks is
constructed using a 1D hydrostatic MARCS model of a RSG star with parameters
matching those of the 3D RHD simulation. The left panel of Fig. 23 compares
the distributions of the line-formation depths (represented by the maxima of
the CF in the $\log\tau_{0}-\rm CF$ plane) in different masks for the 3D
snapshot and the 1D model used for the masks construction. The 1D and 3D line-
formation depths qualitatively overlap in all masks from the inner (mask C1)
to the outer layers (mask C8). Thus, the constructed masks from a 1D model are
usable to probe different atmospheric depths in a dynamical atmosphere.

As the next step, the tomographic masks are cross-correlated with a synthetic
spectrum computed from a 3D snapshot of the same simulation used above. The
resulting cross-correlation functions (CCFs) are shown in the right panel of
Fig. 23 together with the histograms of the line-of-sight ($V_{z}$) velocities
extract from the 3D structure of the snapshot. The velocity field overlaps the
CCF functions across the atmosphere and, as a consequence, the line-of-sight
velocity distribution can be fully recovered. In conclusion, the emerging
spectral lines spread over different optical depths due to the non-radial
convective movements characterizing the stellar atmosphere. This theoretical
validation of the tomographic method opens a new window for the study of
stellar dynamical cycles in evolved stars, and in particular RSGs.

Figure 24: Hysteresis loops between the temperature and RV variations for a
single light cycle of $\mu$ Cep (left, from Kravchenko et al., 2019),
Betelgeuse (middle, from Kravchenko et al., 2021) and the 3D simulation
(right, from Kravchenko et al., 2019), and for the different tomographic masks
(rows). The arrow indicates the direction of evolution along the hysteresis
loops. Figure 25: Velocity maps for the C4 mask and for different 3D snapshots
along the hysteresis loop for that mask ( right panel of Fig. 24). The red and
blue colors correspond to inward and outward moving material, respectively.
The central panel shows the behaviour of the radial velocity. The arrow
indicates the direction of the evolution along the hysteresis loop,
characteristic signature of the convective turn-over of the material in the
stellar atmosphere.

####  4.1.4 Convective cycles in RSG stars: the cases of $\mu$ Cep and
Betelgeuse

The first application of the tomographic method to interpret stellar dynamical
cycles in RSGs is reported by Kravchenko et al. (2019) and Kravchenko et al.
(2021) for two prototypical and extensively observed RSG stars: $\mu$ Cep and
Betelgeuse. In particular, these studies aims at explaining the short-period,
photometric variability of RSGs, whose nature is still debated due to the
semi- and irregular character of the variations.

To identify the acting mechanism(s), the tomographic method is applied to
several years of regular time-series of high-resolution (resolving power of
$\lambda/\Delta\lambda\sim$80000) optical spectra (3800$-$9000 Å) of these
objects, where a phase shift was detected between the radial velocity and
photometric variations (Kravchenko et al., 2019, 2021).
This phase shift results in a hysteresis loop (Gray, 2008) in the temperature-
velocity plane (left and middle panels of Fig. 24) with timescales similar to
the photometric ones.

On the 3D RHD simulation side, its atmosphere reveals qualitatively similar
hysteresis-like behavior in the temperature–velocity plane (right panel of
Fig. 24) with (quasi-) periods consistent with the observed photometric
timescales. In the simulation, the hysteresis process is the characteristic
signature of the convective turn-over of the material in the stellar
atmosphere.
As an example, Figure 25 illustrates the velocity maps for mask C4, for
different 3D snapshots along the hysteresis loop from the right panel of Fig.
24. They reveal upward (blue) and downward (red) motions of matter extending
over large portions of the stellar surface: the relative fraction of upward
and downward motions is what distinguishes the upper from the lower part of
the hysteresis loop, its top part (zero velocity) being characterized by equal
surfaces of rising and falling material. The bottom part of the hysteresis
loop occurs, as expected when the stellar surface is covered mostly by down-
falling material. Moreover, the sound-crossing timescale in the outer layers
is of the same order as the hysteresis-loop timescale. This suggests that
hysteresis loops are linked to acoustic waves originating from the turbulent
convective flow below the visible atmosphere and propagating upward where they
modulate the convective energy flux with shocks.
Another example of the use of this method concerns the RSG Betelgeuse, which
experienced a historic dimming of its visible brightness (whose usual apparent
magnitude comprises between 0.1 and 1.0) down to 1.614$\pm$ 0.008 magnitudes
around 07-13 February 2020 (Guinan et al., 2020). Using this tomography method
and 3D simulations, Kravchenko et al. (2021) showed the presence of two
successive shocks in February 2018 and January 2019 that produced a rapid
expansion of part of the atmosphere of Betelgeuse. The authors provided an
explanation of the phenomenon arguing that the sudden increase in molecular
opacity in the cooler upper atmosphere of Betelgeuse should have caused
obscuration of the star and the observed unusual decrease of its brightness.
Other explanations based on the circumstellar environment where dust grains
form are also possible (e.g. Montargès et al. 2021). For more details, the
reader can refer to the recent review of Wheeler and Chatzopoulos (2023).

###  4.2 Athena++ applications: Calibrating Convection for 1D models and
Predicting Supernova Signatures

The Athena++ simulations have thus far been applied primarily in the context
of 1D stellar and supernova modeling – in particular, stellar convection and
supernova shock breakout. We discuss these applications in the following part
of this Section.

####  4.2.1 Calibrating Mixing Length Theory from the Deep Stellar Interior

Figure 26:  Radial fluid velocity magnitudes at characteristic snapshots of
two Athena++ RSG models compared to MLT (RSG1L4.5, day 4707, left; RSG2L4.9,
day 4927, right). Mass-weighted average velocities are shown as thick black
lines, with 68% and 80% of the mass lying within the dark and light shaded
regions, respectively. Greyscale lines indicate convective velocities
predicted from MLT given the volume-averaged temperature and density profiles
and the model luminosity, for integer values of $\alpha=$1—10 where higher
values of $\alpha$ predict larger convective velocities. The plots are
truncated where the turbulent motions do not resemble MLT-like convection,
which lies inside the region where $\nabla\langle s\rangle\leq 0$. Figure from
Goldberg et al. (2022a).

Convection is a fundamentally 3D process, yet it is a dominant mode of energy
transport and compositional mixing in stellar interiors and envelopes,
necessitating prescriptions for its effects on 1D stellar evolution models.
The Mixing Length Theory (MLT), first formulated for stellar interiors by
Böhm-Vitense (1958) and expanded by Henyey et al. (1965) to model the
transition from convective to radiative transport, serves as the default
prescription for convection in the vast majority of 1D stellar evolution codes
(see Joyce and Tayar, 2023, for a detailed review).

A key element of the MLT prescription is that the characteristic length scale
of transport which determines the efficiency of convection, i.e. the mixing
length $\ell$, taken locally to be a free parameter $\alpha$ times the
pressure scale height, $\ell=\alpha H$. Especially if the _envelope_ is
convective, assumptions about $\alpha$ and other MLT parameters strongly
influence the stellar radii and $T_{\rm eff}$ output by 1D codes (e.g.
Stothers and Chin, 1995; Massey and Olsen, 2003; Meynet and Maeder, 1997;
Meynet et al., 2015; Joyce and Tayar, 2023) (in addition to assumptions
relevant to the structure and location of internal convective boundaries).
Such values of $\alpha$ are often observationally calibrated by comparing
stellar models to the Sun (e.g. Charbonnel and Lebreton, 1993; Cinquegrana and
Joyce, 2022, and many others) and other well-characterized stars (e.g.
Guenther and Demarque, 2000; Chun et al., 2018; Joyce and Chaboyer, 2018).

Computational RHD models of convection enable a theoretical calibration of MLT
free parameters. Early 2D RHD simulations (Ludwig et al., 1999) and 3D
simulations by Sonoi et al. (2019) calibrated MLT in $L\ll L_{\mathrm{Edd}}$
stars, constraining $\alpha\approx 1-2$ across a grid of low-luminosity cool
giant atmospheres (Red Giants with $T_{\rm eff}>4000$K), in agreement with the
observational constraints. Other works (e.g. Trampedach et al. 2014; Magic et
al. 2015b; Salaris and Cassisi 2015) recovered similar calibrations in similar
$L\ll L_{\mathrm{Edd}}$ stars, with some trends towards higher efficiency
(i.e. higher $\alpha$) in more luminous stars. As convection becomes even more
vigorous in the coolest, most luminous stars, such as AGB and RSG stars,
individual convective plumes take up larger and larger fractions of the star
and global or semi-global simulations in full-star geometry are required, and
in many cases convective motions particularly at the stellar surface begin to
deviate significantly from MLT (see, e.g., the discussion in Trampedach et
al., 2013).

Figure 27:  Left: Comparison of superadiabaticity (upper left) and specific
entropy (lower left) derived from the RSG2L4.9 Athena++ RSG model (black
dashed lines) and from MLT (solid colored lines) when $P_{\rm turb}$ is
neglected (i.e. when $P=P_{\rm therm}$ and
$\nabla\equiv\nabla_{\mathrm{therm}}=d\ln T/d\ln P_{\rm therm}$). Right:
Impact of turbulent pressure on the mixing length calibration. Comparison of
superadiabaticity (middle right) and entropy (lower right) are shown, for our
averaged 3D RSG2L4.9 model (black dashed lines) and MLT with different
$\alpha$ (solid colored lines) as a function of temperature. The upper right
panel shows $P_{\rm turb}/P_{\rm therm}$ (dotted line) and $d\ln P_{\rm
therm}/d\ln P_{\rm tot}$ (dash-dot line), which following Henyey et al. (1965)
are then included in the MLT calculation and in the averaged 3D $\nabla$. The
grey dashed line in the middle right panel gives $d\ln T/d\ln P_{\rm therm}$,
matching the black dashed line on the upper left panel. All values shown are
derived from the time-averaged, shell-averaged density and temperature
profiles, as well as the time-averaged luminosity at the simulation outer
boundary, beyond day 4500 in RSG2L4.9. Figures adapted from Goldberg et al.
(2022a).

Using the Athena++ setup described above, Goldberg et al. (2022a) produced the
first 3D RHD calibration of MLT in the RSG regime. Despite large deviations
from MLT near the surface, there exists a large region in the Athena++ RSG
envelopes where the convective motions and fluxes do resemble those described
by MLT (i.e. where convection carries a majority of the flux, where the radial
entropy gradient is zero or negative, and where positive density fluctuations
correlate with inwards radial velocities). Therefore the convective velocities
and thermodynamic gradients (in particular $\nabla\equiv d\ln T/d\ln P$ and
the entropy) which are recovered self-consistently from the 3D RHD equations
can be directly compared against analytic expectations from MLT in order to
calibrate $\alpha$. The velocity profiles were found to be flatter than MLT
would predict, likely owing to the non-local nature of the coherent convective
plumes, with good order-of-magnitude agreement but large scatter exceeding a
factor of 10 variation of the mixing length parameter $\alpha$ (Fig. 26).
Comparing $\nabla$ and the entropy with different $\alpha$, Goldberg et al.
(2022a) recovered high convective efficiencies ($\alpha\approx 3-4$). This
comparison is shown in Fig. 27 for calibrations both without (left) and
including (right) turbulent pressure via the method suggested by Henyey et al.
(1965) but oft-neglected in 1D stellar evolution codes. This is consistent
with a number of estimates of larger-than-solar values of $\alpha$ in the RSG
regime, from observed RSG populations on the HR diagram compared to stellar
models (Chun et al., 2018) and early supernova colors (Dessart et al., 2013).

####  4.2.2 Modeling Shock Breakout Emission in Type II-P Supernovae

Another novel application of the Athena++ framework is the explosion of these
large-scale convective 3D models. As the real explosion would occur beneath
the inner boundary of the Athena++ simulation domain, in the setup detailed in
Goldberg et al. (2022b), an initially spherical shock is driven through the
Athena++ simulation inner boundary, fed by a 1D thermal bomb explosion deep in
a MESA stellar model.

Figure 28: 3D surface rendering of the radius where $\tau_{s}=100$ along each
radial direction for the RSG1L4.5 model, where SBO is expected to occur. Four
characteristic lines of sight (A,B,C,D) are labeled corresponding to a
“valley”, “hillside”, “plateau”, and “mountain”, of the SBO surface.

In the death of a RSG as a Type II (Hydrogen-rich) Supernova, a strong
radiation-dominated shock propagates through the convective envelope. As the
shock approaches the outer layers, it reaches optical depth sufficiently low
for photons carrying the explosive luminosity to escape ahead of the emerging
shock via diffusion, at $\tau\lesssim c/v_{\rm sh}$ where $c$ is the speed of
light and $v_{\rm sh}$ is the shock velocity. This leads to a hot
($T>10^{5}$K), bright flash, known as the “shock breakout” (SBO). There is an
abundance of predictions for the brightness and duration of SBO in spherical-
symmetry, (e.g. Nakar and Sari 2010; Rabinak and Waxman 2011; Sapir et al.
2011; Katz et al. 2012; Sapir et al. 2013; Shussman et al. 2016; Sapir and
Waxman 2017; Kozyreva et al. 2020; Morag et al. 2023). One prediction of the
1D semi-analytics is that the observed SBO duration may be set by the light-
travel time over the stellar radius, $R/c\lesssim 1$ hour ($\approx$30 min for
$800R_{\odot}$), which if measured would provide a direct constraint on the
stellar radius. This timescale has not yet been seen in observations, which
are currently limited to a few serendipitous detections from the GALEX
mission, all of which show durations of $\gtrsim 6$ hours (e.g. Schawinski et
al., 2008; Gezari et al., 2008, 2010, 2015). The data are expected to improve
significantly with future satellite missions such as ULTRASAT (Sagiv et al.,
2014; Asif et al., 2021; Shvartzvald et al., 2023), which is expected to
capture hundreds of SBO’s at high cadence.

The 3D envelopes differ from standard 1D models in two important ways. First
is the extended halo of low-density material outside the traditional stellar
photosphere, which becomes ionized and contributes non-negligibly to the
optical depth as the supernova shock approaches. The turbulent pressure in
this region extends the scale height of the low-density outer layers discussed
in the context of 1D models (e.g. Morag et al., 2023), resembling the
exceedingly low densities required to reproduce SBO observations (Schawinski
et al., 2008). The second is the large-scale intrinsically 3D convective
granulation, spanning large fractions of the stellar surface and spanning
“valley-to-mountain depths” of $\Delta R\approx 10-20\%$ of the (large)
stellar radius in both the CO5BOLD and Athena++ simulations, shown in Figure
28. This global structure is also consistent with tomographic observations as
discussed in Section 4.1.

Figure 29: Gas temperature of the fastest-moving ejecta along each radial line
of sight (left column) and the radiative luminosity $4\pi r^{2}F_{r}$ at the
approximate photospheric location along each radial ray for a given
($\theta,\phi$), defined as the first location where the radiation
flux/radiation energy density $F_{r}/E_{r}=c/3$ inwards along each radial ray.
The temperature of the fastest-moving ejecta signals whether thermal energy is
trapped behind the shock pre-“breakout” (indicated by hotter temperatures) or
whether it is able to escape and cool (cooler temperatures). Time until peak
bolometric luminosity is shown to the right of each row, and the panels are
sampled approximately one diffusion-time apart. Reference locations A, B, C,
and D indicated in Fig. 28 are shown in the ($T_{\mathrm{fast}}$, -3.5 hours)
panel. Figure adapted from Goldberg et al. (2022b). Figure 30: Bolometric
lightcurves ($L_{\mathrm{bol}}$) calculated for varied progenitor model (top)
and explosion energies (bottom). Black curves correspond to an $8\times
10^{50}$erg explosion of the RSG1L4.5 model (with $\Delta R=80R_{\odot}$).
Faint lines (calculated at $r=2700R_{\odot}$) show the integrated flux an
external observer would see looking back at the star from a particular
$\theta,\phi$ location (note that this is not shown at the photosphere as in
Fig. 29, but rather for an external observer viewing the whole explosion at
different angles from far away). For comparison, SBO lightcurves from
spherically symmetric MESA+STELLA explosion models (Goldberg and Bildsten,
2020) are shown in the upper panel (blues), labeled by ejecta
mass/$M_{\odot}$, radius/$R_{\odot}$, and explosion
energy/$10^{51}\mathrm{erg}$, where the model truncates at the stellar
photosphere. The x-axis for all 3D curves is the time since the maximum shell-
averaged $L_{\rm bol}$ (i.e. the peak of the thick curves). Figure from
Goldberg et al. (2022b).

Goldberg et al. (2022b) find that there are three timescales which can
influence SBO. The first timescale is the diffusion time $t_{\rm diff}$ for
photons to leak out from the $\tau=c/v_{\rm sh}$ location (e.g. Nakar and
Sari, 2010), which is somewhat prolonged by the presence of the low-density
halo above the photosphere. The second timescale is $R/c$, which will smear
out the SBO signal if larger than $t_{\rm diff}$ (e.g. Katz et al., 2012). The
third timescale is introduced by the large-scale 3D granulation which shapes
the surface where shock breakout occurs. This can be seen in Figure 28 showing
radii of constant $\tau_{s}=100$ corresponding roughly to the breakout surface
$c/v_{\rm sh}$ for shock velocities of a few thousand km/s, where
$\tau_{s}=\int_{r}^{\infty}\kappa_{s}\rho dr$ is the optical depth to electron
scattering (with opacity $\kappa_{s}\approx 0.32$ cm2 g-1). Four points of
varying topography on this surface are labeled A,B,C,D ascending from deepest
to shallowest. For the RSG1L4.5 envelope shown, the large-scale inhomogeneity
entails radial differences in the SBO surface topography spanning $\Delta
R\approx 80R_{\odot}$ (for RSG2L4.9, $\Delta R\approx 200R_{\odot}$). As the
shock traverses this large $\Delta R$, it breaks out at different radii at
different times along different patches of the simulation domain, first in the
“valleys” (e.g. A) and lastly along the “mountains” (e.g. D). The third
timescale is thus the time for the shock to cross $\Delta R$: $t_{\rm
cross}\approx\Delta R/v_{\rm sh}$. When $v_{\rm sh}$ is sufficiently large,
$v_{\rm sh}>c/(\kappa\rho_{0}\Delta R)\approx 1300$ km/s for $\Delta R\approx
80R_{\odot}$ and a typical breakout density $\rho_{0}\approx 1.25\times
10^{-10}\,\mathrm{g\,cm^{-3}}$ and $\kappa\approx
0.32\mathrm{cm}^{2}\,\mathrm{g}^{-1}$, SBO from different patches on the
stellar surface at different times will smear out the observed SBO duration.
This effect is seen in the evolution of the ejecta temperature and radiative
flux at a given patch shown in (Figure 29); “valleys” (e.g. A) undergo SBO
earliest, and even at peak L, some “mountains” (e.g. D) are still pre-SBO,
still hot and comparatively dark. Lightcurves (Figure 30) from 3D envelope
explosions show a slower rise to peak brightness, with a gradual rise even
before half-peak-brightness occurring over multiple hours, with longer
durations from half-peak-brightness to half-peak-brightness spanning $\Delta
t_{1/2}\approx 3-6$ hours. Correspondingly, the predicted luminosities are
lowered by a factor of $2-10$ (to $L_{\mathrm{bol}}\sim 10^{44}\mathrm{erg\
s^{-1}}$) as comparable radiation energy escapes over a longer timescale.

Because $t_{\rm cross}$ exceeds $t_{\rm diff}$ except in the limit of low
shock velocities and thereby large diffusion times for large $\Delta R$
($\gtrsim 10$% of the stellar radius in both CO5BOLD and Athena++
simulations), Goldberg et al. (2022b) argue that in either case, $t_{\rm
diff}$ and $t_{\rm cross}$ will both exceed $R/c$, rendering it unlikely to
measure RSG radii from the SBO duration. However, this introduces the
possibility that future observations of SBO can be used to constrain $\Delta
R$ and provide yet another probe into the structure and inhomogeneity of the
outer layers of cool evolved stars.

##  5 Summary

Convection is a non-local and complex process with non-linear interactions
over many disparate length scales happening in a 3-dimensional geometry. To
have a full understanding of this phenomenon, global 3D radiation-
hydrodynamics simulations are crucial to modeling the convective atmospheres
of RSG and AGB stars. These simulations avoid the use of parameterized
convection (e.g., in terms of mixing length theory, commonly used in classic,
1D, linear and non-linear pulsation models) and make it possible to follow the
flowing matter, in complete 3D geometry, from the stellar interior to the
circumstellar envelope, via the stellar atmosphere and the dust forming
region.

3D RHD simulations predict how long-lasting giant convection cells (with sizes
comparable to the stellar radii), consisting of wide upflow regions that are
surrounded by turbulent downdrafts, together with short-lived smaller surface
granules and radial fundamental-mode pulsations, give rise to variable surface
structures that, as a consequence, affect the triggering and shape of the
strong stellar winds as well as their fundamental parameters. These 3D
insights can be contextualized and aid in the development and calibration of
further 1D models.

From an observational point of view, the dynamical processes across the
convective atmospheres largely impact the formation of spectral lines with
consequences on photometry, radial velocities, abundances. The associated
intrinsic and physical variability can be used to recover the projected
velocity field (along the line-of-sight) at different atmospheric layers using
the tomographic method and 3D simulations. This theoretical interpretation
suggests that the detected hysteresis loops are the characteristic signal of
stellar convection and are intimately linked with the acoustic waves
originating from the turbulent convective flow propagating upward to the
surface layers where the convective energy is modulated by supersonic shocks.
The application of tomography with 3D simulations is a powerful technique for
a comprehensive study of the stellar dynamical cycles in evolved stars.

###### Acknowledgements.

The authors thank the referees for the detailed and constructive comments to
improve this work. This work is supported by the French National Research
Agency (ANR) funded project PEPPER (ANR-20-CE31-0002). This work was granted
access to the HPC resources of Observatoire de la Côte d’Azur $--$ Mésocentre
SIGAMM and resources provided by the Swedish National Infrastructure for
Computing (SNIC) at UPPMAX. HPC resources supported by the NASA ATP grant
ATP-80NSSC18K0560 were provided by the NASA High-End Computing (HEC) program
through the NASA Advanced Supercomputing (NAS) Division at Ames Research
Center. The Flatiron Institute is supported by the Simons Foundation.

## References

  * Ahmad et al. (2023) Ahmad A, Freytag B, Höfner S (2023) Properties of self-excited pulsations in 3D simulations of AGB stars and red supergiants. A&A669:A49. https://doi.org/10.1051/0004-6361/202244555. arXiv:2211.07682 [astro-ph.SR]
  * Albrow and Cottrell (1996) Albrow MD, Cottrell PL (1996) Contribution functions and the depths of formation of spectral lines in Cepheids. MNRAS278:337. https://doi.org/10.1093/mnras/278.2.337
  * Allende Prieto et al. (2013) Allende Prieto C, Koesterke L, Ludwig HG, Freytag B, Caffau E (2013) Convective line shifts for the Gaia RVS from the CIFIST 3D model atmosphere grid. A&A 550:A103. https://doi.org/10.1051/0004-6361/201220064. arXiv:1301.3703 [astro-ph.SR]
  * Alvarez et al. (2000) Alvarez R, Jorissen A, Plez B, Gillet D, Fokin A (2000) Envelope tomography of long-period variable stars. I. The Schwarzschild mechanism and the Balmer emission lines. A&A362:655–665. https://doi.org/10.48550/arXiv.astro-ph/0010132. arXiv:astro-ph/0010132 [astro-ph]
  * Alvarez et al. (2001) Alvarez R, Jorissen A, Plez B, et al. (2001) Envelope tomography of long-period variable stars II. Method. A&A379:288–304. https://doi.org/10.1051/0004-6361:20011294. arXiv:astro-ph/0109251 [astro-ph]
  * Antoni and Quataert (2022) Antoni A, Quataert E (2022) Numerical simulations of the random angular momentum in convection: Implications for supergiant collapse to form black holes. MNRAS511(1):176–197. https://doi.org/10.1093/mnras/stab3776. arXiv:2107.09068 [astro-ph.HE]
  * Asif et al. (2021) Asif A, Barschke M, Bastian-Querner B, et al. (2021) Design of the ULTRASAT UV camera. In: Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series. Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, vol 11821. p 118210U. https://doi.org/10.1117/12.2594253. arXiv:2108.01493 [astro-ph.IM]
  * Asplund et al. (2000) Asplund M, Ludwig HG, Nordlund Å, Stein RF (2000) The effects of numerical resolution on hydrodynamical surface convection simulations and spectral line formation. A&A 359:669–681. arXiv:astro-ph/0005319 [astro-ph]
  * Asplund et al. (2009) Asplund M, Grevesse N, Sauval AJ, Scott P (2009) The Chemical Composition of the Sun. A&AR 47:481–522. https://doi.org/10.1146/annurev.astro.46.060407.145222. arXiv:0909.0948
  * Aurière et al. (2010) Aurière M, Donati JF, Konstantinova-Antova R, et al. (2010) The magnetic field of Betelgeuse: a local dynamo from giant convection cells? A&A 516:L2. https://doi.org/10.1051/0004-6361/201014925. arXiv:1005.4845 [astro-ph.SR]
  * Bigot et al. (2011) Bigot L, Mourard D, Berio P, et al. (2011) The diameter of the CoRoT target HD 49933. Combining the 3D limb darkening, asteroseismology, and interferometry. A&A534:L3. https://doi.org/10.1051/0004-6361/201117349. arXiv:1110.0985 [astro-ph.SR]
  * Böhm-Vitense (1958) Böhm-Vitense E (1958) Über die Wasserstoffkonvektionszone in Sternen verschiedener Effektivtemperaturen und Leuchtkräfte. Mit 5 Textabbildungen. ZAp46:108
  * Bonifacio et al. (2018) Bonifacio P, Caffau E, Ludwig HG, et al. (2018) Using the CIFIST grid of CO5BOLD 3D model atmospheres to study the effects of stellar granulation on photometric colours. I. Grids of 3D corrections in the UBVRI, 2MASS, HIPPARCOS, Gaia, and SDSS systems. A&A611:A68. https://doi.org/10.1051/0004-6361/201732232. arXiv:1712.00024 [astro-ph.SR]
  * Caffau et al. (2011) Caffau E, Ludwig HG, Steffen M, Freytag B, Bonifacio P (2011) Solar Chemical Abundances Determined with a CO5BOLD 3D Model Atmosphere. Sol. Phys.268(2):255–269. https://doi.org/10.1007/s11207-010-9541-4. arXiv:1003.1190 [astro-ph.SR]
  * Cegla (2019) Cegla H (2019) The Impact of Stellar Surface Magnetoconvection and Oscillations on the Detection of Temperate, Earth-Mass Planets Around Sun-Like Stars. Geosciences 9(3):114. https://doi.org/10.3390/geosciences9030114. arXiv:1904.03200 [astro-ph.EP]
  * Charbonnel and Lebreton (1993) Charbonnel C, Lebreton Y (1993) Standard solar model: Interplay between the equation of state, the opacity and the determination of the initial helium content. A&A280(2):666–674
  * Chiavassa and Brogi (2019) Chiavassa A, Brogi M (2019) Planet and star synergy at high-spectral resolution. A rationale for the characterization of exoplanet atmospheres. I. The infrared. A&A631:A100. https://doi.org/10.1051/0004-6361/201936566. arXiv:1909.11807 [astro-ph.EP]
  * Chiavassa and Freytag (2015) Chiavassa A, Freytag B (2015) 3D Hydrodynamical Simulations of Evolved Stars and Observations of Stellar Surfaces. In: Kerschbaum F, Wing RF, Hron J (eds) Why Galaxies Care about AGB Stars III: A Closer Look in Space and Time. Astronomical Society of the Pacific Conference Series, vol 497. p 11. arXiv:1410.3868 [astro-ph.SR]
  * Chiavassa et al. (2009) Chiavassa A, Plez B, Josselin E, Freytag B (2009) Radiative hydrodynamics simulations of red supergiant stars. I. interpretation of interferometric observations. A&A 506(3):1351–1365. https://doi.org/10.1051/0004-6361/200911780. arXiv:0907.1860 [astro-ph.SR]
  * Chiavassa et al. (2010) Chiavassa A, Haubois X, Young JS, et al. (2010) Radiative hydrodynamics simulations of red supergiant stars. II. Simulations of convection on Betelgeuse match interferometric observations. A&A 515:A12. https://doi.org/10.1051/0004-6361/200913907. arXiv:1003.1407 [astro-ph.SR]
  * Chiavassa et al. (2011a) Chiavassa A, Freytag B, Masseron T, Plez B (2011a) Radiative hydrodynamics simulations of red supergiant stars. IV. Gray versus non-gray opacities. A&A 535:A22. https://doi.org/10.1051/0004-6361/201117463. arXiv:1109.3619 [astro-ph.SR]
  * Chiavassa et al. (2011b) Chiavassa A, Pasquato E, Jorissen A, et al. (2011b) Radiative hydrodynamic simulations of red supergiant stars. III. Spectro-photocentric variability, photometric variability, and consequences on Gaia measurements. A&A 528:A120. https://doi.org/10.1051/0004-6361/201015768. arXiv:1012.5234 [astro-ph.SR]
  * Chiavassa et al. (2018a) Chiavassa A, Casagrande L, Collet R, et al. (2018a) The STAGGER-grid: A grid of 3D stellar atmosphere models. V. Synthetic stellar spectra and broad-band photometry. A&A 611:A11. https://doi.org/10.1051/0004-6361/201732147. arXiv:1801.01895 [astro-ph.SR]
  * Chiavassa et al. (2018b) Chiavassa A, Freytag B, Schultheis M (2018b) Heading Gaia to measure atmospheric dynamics in AGB stars. A&A 617:L1. https://doi.org/10.1051/0004-6361/201833844. arXiv:1808.02548 [astro-ph.SR]
  * Chiavassa et al. (2020) Chiavassa A, Kravchenko K, Millour F, et al. (2020) Optical interferometry and Gaia measurement uncertainties reveal the physics of asymptotic giant branch stars. A&A 640:A23. https://doi.org/%**** ̵evolved˙stars˙convection˙review˙v5.bbl ̵Line ̵200 ̵****10.1051/0004-6361/202037832. arXiv:2006.07318 [astro-ph.SR]
  * Chiavassa et al. (2021) Chiavassa A, Kravchenko K, Montargès M, et al. (2021) The extended atmosphere and circumstellar environment of the cool evolved star VX Sagittarii as seen by MATISSE. arXiv e-prints arXiv:2112.10695. arXiv:2112.10695 [astro-ph.SR]
  * Chiavassa et al. (2022a) Chiavassa A, Kravchenko K, Montargès M, et al. (2022a) The extended atmosphere and circumstellar environment of the cool evolved star VX Sagittarii as seen by MATISSE. A&A658:A185. https://doi.org/10.1051/0004-6361/202142514. arXiv:2112.10695 [astro-ph.SR]
  * Chiavassa et al. (2022b) Chiavassa A, Kudritzki R, Davies B, Freytag B, de Mink SE (2022b) Probing red supergiant dynamics through photo-center displacements measured by Gaia. A&A661:L1. https://doi.org/10.1051/0004-6361/202243568. arXiv:2205.05156 [astro-ph.SR]
  * Chun et al. (2018) Chun SH, Yoon SC, Jung MK, Kim DU, Kim J (2018) Evolutionary Models of Red Supergiants: Evidence for A Metallicity-dependent Mixing Length and Implications for Type IIP Supernova Progenitors. ApJ853(1):79. https://doi.org/10.3847/1538-4357/aa9a37. arXiv:1711.03965 [astro-ph.SR]
  * Cinquegrana and Joyce (2022) Cinquegrana GC, Joyce M (2022) Solar Calibration of the Convective Mixing Length for Use with the ÆSOPUS Opacities in MESA. Research Notes of the American Astronomical Society 6(4):77. https://doi.org/10.3847/2515-5172/ac6611
  * Cranmer and Saar (2011) Cranmer SR, Saar SH (2011) Testing a Predictive Theoretical Model for the Mass Loss Rates of Cool Stars. ApJ 741(1):54. https://doi.org/10.1088/0004-637X/741/1/54. arXiv:1108.4369 [astro-ph.SR]
  * Creevey et al. (2012) Creevey OL, Thévenin F, Boyajian TS, et al. (2012) Fundamental properties of the Population II fiducial stars ¡ASTROBJ¿HD 122563¡/ASTROBJ¿ and ¡ASTROBJ¿Gmb 1830¡/ASTROBJ¿ from CHARA interferometric observations. A&A545:A17. https://doi.org/10.1051/0004-6361/201219651. arXiv:1207.5954 [astro-ph.SR]
  * De Beck et al. (2010) De Beck E, Decin L, de Koter A, et al. (2010) Probing the mass-loss history of AGB and red supergiant stars from CO rotational line profiles. II. CO line survey of evolved stars: derivation of mass-loss rate formulae. A&A 523:A18. https://doi.org/10.1051/0004-6361/200913771. arXiv:1008.1083 [astro-ph.SR]
  * Decin et al. (2020) Decin L, Montargès M, Richards AMS, et al. (2020) (Sub)stellar companions shape the winds of evolved stars. Science 369(6510):1497–1500. https://doi.org/10.1126/science.abb1229. arXiv:2009.11694 [astro-ph.SR]
  * Dessart et al. (2013) Dessart L, Hillier DJ, Waldman R, Livne E (2013) Type II-Plateau supernova radiation: dependences on progenitor and explosion properties. MNRAS433:1745–1763. https://doi.org/10.1093/mnras/stt861. arXiv:1305.3386 [astro-ph.SR]
  * Freytag (2013) Freytag B (2013) Advances in the hydrodynamics solver of CO5BOLD. Memorie della Societa Astronomica Italiana Supplementi 24:26
  * Freytag (2017) Freytag B (2017) Boundary conditions in CO5BOLD. Mem. Soc. Astron. Italiana88:12
  * Freytag and Höfner (2023) Freytag B, Höfner S (2023) Global 3D radiation-hydrodynamical models of AGB stars with dust-driven winds. A&A669:A155. https://doi.org/10.1051/0004-6361/202244992. arXiv:2301.11836 [astro-ph.SR]
  * Freytag et al. (2002) Freytag B, Steffen M, Dorch B (2002) Spots on the surface of Betelgeuse – Results from new 3D stellar convection models. Astronomische Nachrichten 323:213–219. https://doi.org/10.1002/1521-3994(200208)323:3/4¡213::AID-ASNA213¿3.0.CO;2-H
  * Freytag et al. (2012) Freytag B, Steffen M, Ludwig HG, et al. (2012) Simulations of stellar convection with CO5BOLD. Journal of Computational Physics 231(3):919–959. https://doi.org/10.1016/j.jcp.2011.09.026. arXiv:1110.6844 [astro-ph.SR]
  * Freytag et al. (2017) Freytag B, Liljegren S, Höfner S (2017) Global 3D radiation-hydrodynamics models of AGB stars. Effects of convection and radial pulsations on atmospheric structures. A&A 600:A137. https://doi.org/10.1051/0004-6361/201629594. arXiv:1702.05433 [astro-ph.SR]
  * Gezari et al. (2008) Gezari S, Dessart L, Basa S, et al. (2008) Probing Shock Breakout with Serendipitous GALEX Detections of Two SNLS Type II-P Supernovae. ApJ683(2):L131. https://doi.org/10.1086/591647. arXiv:0804.1123 [astro-ph]
  * Gezari et al. (2010) Gezari S, Rest A, Huber ME, et al. (2010) GALEX and Pan-STARRS1 Discovery of SN IIP 2010aq: The First Few Days After Shock Breakout in a Red Supergiant Star. ApJ720(1):L77–L81. https://doi.org/10.1088/2041-8205/720/1/L77. arXiv:1007.4551 [astro-ph.CO]
  * Gezari et al. (2015) Gezari S, Jones DO, Sanders NE, et al. (2015) GALEX Detection of Shock Breakout in Type IIP Supernova PS1-13arp: Implications for the Progenitor Star Wind. ApJ804(1):28. https://doi.org/10.1088/0004-637X/804/1/28. arXiv:1502.06964 [astro-ph.HE]
  * Goldberg and Bildsten (2020) Goldberg JA, Bildsten L (2020) The Value of Progenitor Radius Measurements for Explosion Modeling of Type II-Plateau Supernovae. ApJ895(2):L45. https://doi.org/10.3847/2041-8213/ab9300. arXiv:2005.07290 [astro-ph.SR]
  * Goldberg et al. (2022a) Goldberg JA, Jiang YF, Bildsten L (2022a) Numerical simulations of convective three-dimensional red supergiant envelopes. ApJ929(2):156. https://doi.org/10.3847/1538-4357/ac5ab3, URL https://doi.org/10.3847/1538-4357/ac5ab3
  * Goldberg et al. (2022b) Goldberg JA, Jiang YF, Bildsten L (2022b) Shock Breakout in Three-dimensional Red Supergiant Envelopes. ApJ933(2):164. https://doi.org/10.3847/1538-4357/ac75e3. arXiv:2206.04134 [astro-ph.SR]
  * Gray (2008) Gray DF (2008) Mass Motions in the Photosphere of Betelgeuse. AJ135(4):1450–1458. https://doi.org/10.1088/0004-6256/135/4/1450
  * Guenther and Demarque (2000) Guenther DB, Demarque P (2000) $\alpha$ Centauri AB. ApJ531(1):503–520. https://doi.org/10.1086/308446
  * Guinan et al. (2020) Guinan E, Wasatonic R, Calderwood T, Carona D (2020) The Fall and Rise in Brightness of Betelgeuse. The Astronomer’s Telegram 13512:1
  * Gustafsson et al. (2008) Gustafsson B, Edvardsson B, Eriksson K, et al. (2008) A grid of MARCS model atmospheres for late-type stars. I. Methods and general properties. A&A486(3):951–970. https://doi.org/10.1051/0004-6361:200809724. arXiv:0805.0554 [astro-ph]
  * Harten et al. (1983) Harten A, Lax P, Leer B (1983) Development and Application of Two-Dimensional Numerical Model on Shallow Water Flows Using Finite Volume Method. SIAM Review 25:35
  * Haubois et al. (2009) Haubois X, Perrin G, Lacour S, et al. (2009) Imaging the spotty surface of Betelgeuse in the H band. A&A 508(2):923–932. https://doi.org/10.1051/0004-6361/200912927. arXiv:0910.4167 [astro-ph.SR]
  * Haubois et al. (2019) Haubois X, Norris B, Tuthill PG, et al. (2019) The inner dust shell of Betelgeuse detected by polarimetric aperture-masking interferometry. A&A628:A101. https://doi.org/10.1051/0004-6361/201833258. arXiv:1907.08594 [astro-ph.SR]
  * Hauschildt et al. (1997) Hauschildt PH, Baron E, Allard F (1997) Parallel Implementation of the PHOENIX Generalized Stellar Atmosphere Program. ApJ483(1):390–398. https://doi.org/10.1086/304233. arXiv:astro-ph/9607087 [astro-ph]
  * Henyey et al. (1965) Henyey L, Vardya MS, Bodenheimer P (1965) Studies in Stellar Evolution. III. The Calculation of Model Envelopes. ApJ142:841. https://doi.org/10.1086/148357
  * Höfner and Freytag (2019) Höfner S, Freytag B (2019) Exploring the origin of clumpy dust clouds around cool giants. A global 3D RHD model of a dust-forming M-type AGB star. A&A623:A158. https://doi.org/10.1051/0004-6361/201834799. arXiv:1902.04074 [astro-ph.SR]
  * Höfner and Olofsson (2018) Höfner S, Olofsson H (2018) Mass loss of stars on the asymptotic giant branch. Mechanisms, models and measurements. A&AR 26(1):1. https://doi.org/10.1007/s00159-017-0106-5
  * Iglesias and Rogers (1996) Iglesias CA, Rogers FJ (1996) Updated Opal Opacities. ApJ464:943. https://doi.org/10.1086/177381
  * Iglesias et al. (1992) Iglesias CA, Rogers FJ, Wilson BG (1992) Spin-Orbit Interaction Effects on the Rosseland Mean Opacity. ApJ397:717. https://doi.org/10.1086/171827
  * Jermyn et al. (2022) Jermyn AS, Anders EH, Lecoanet D, Cantiello M, Goldberg JA (2022) Measures of Efficiency of Convection. Research Notes of the American Astronomical Society 6(2):29. https://doi.org/10.3847/2515-5172/ac531e
  * Jermyn et al. (2023) Jermyn AS, Bauer EB, Schwab J, et al. (2023) Modules for Experiments in Stellar Astrophysics (MESA): Time-dependent Convection, Energy Conservation, Automatic Differentiation, and Infrastructure. ApJS265(1):15. https://doi.org/10.3847/1538-4365/acae8d. arXiv:2208.03651 [astro-ph.SR]
  * Jiang (2021) Jiang YF (2021) An Implicit Finite Volume Scheme to Solve the Time-dependent Radiation Transport Equation Based on Discrete Ordinates. ApJS253(2):49. https://doi.org/%**** ̵evolved˙stars˙convection˙review˙v5.bbl ̵Line ̵475 ̵****10.3847/1538-4365/abe303. arXiv:2102.02212 [astro-ph.IM]
  * Jiang (2022) Jiang YF (2022) Multigroup Radiation Magnetohydrodynamics Based on Discrete Ordinates including Compton Scattering. ApJS263(1):4. https://doi.org/10.3847/1538-4365/ac9231. arXiv:2209.06240 [astro-ph.IM]
  * Jiang et al. (2014) Jiang YF, Stone JM, Davis SW (2014) An Algorithm for Radiation Magnetohydrodynamics Based on Solving the Time-dependent Transfer Equation. ApJS213(1):7. https://doi.org/10.1088/0067-0049/213/1/7. arXiv:1403.6126 [astro-ph.IM]
  * Jiang et al. (2015) Jiang YF, Cantiello M, Bildsten L, Quataert E, Blaes O (2015) Local Radiation Hydrodynamic Simulations of Massive Star Envelopes at the Iron Opacity Peak. ApJ813(1):74. https://doi.org/10.1088/0004-637X/813/1/74. arXiv:1509.05417 [astro-ph.SR]
  * Jiang et al. (2018) Jiang YF, Cantiello M, Bildsten L, et al. (2018) Outbursts of luminous blue variable stars from variations in the helium opacity. Nature561(7724):498–501. https://doi.org/10.1038/s41586-018-0525-0. arXiv:1809.10187 [astro-ph.SR]
  * Jorissen et al. (2016) Jorissen A, Van Eck S, Kravchenko K (2016) Velocity Fields in Stellar Atmospheres Probed by Tomography. In: Boffin HMJ, Hussain G, Berger JP, Schmidtobreick L (eds) Astronomy at High Angular Resolution. Astrophysics and Space Science Library, vol 439. p 137. https://doi.org/10.1007/978-3-319-39739-9˙8
  * Josselin and Plez (2007) Josselin E, Plez B (2007) Atmospheric dynamics and the mass loss process in red supergiant stars. A&A469(2):671–680. https://doi.org/10.1051/0004-6361:20066353. arXiv:0705.0266 [astro-ph]
  * Joyce and Chaboyer (2018) Joyce M, Chaboyer B (2018) Not All Stars Are the Sun: Empirical Calibration of the Mixing Length for Metal-poor Stars Using One-dimensional Stellar Evolution Models. ApJ856(1):10. https://doi.org/10.3847/1538-4357/aab200. arXiv:1712.05082 [astro-ph.SR]
  * Joyce and Tayar (2023) Joyce M, Tayar J (2023) A Review of the Mixing Length Theory of Convection in 1D Stellar Modeling. Galaxies 11(3):75. https://doi.org/10.3390/galaxies11030075. arXiv:2303.09596 [astro-ph.SR]
  * Katz et al. (2012) Katz B, Sapir N, Waxman E (2012) Non-relativistic Radiation Mediated Shock Breakouts. II. Bolometric Properties of Supernova Shock Breakout. ApJ747(2):147. https://doi.org/10.1088/0004-637X/747/2/147. arXiv:1103.5276 [astro-ph.HE]
  * Kee et al. (2021) Kee ND, Sundqvist JO, Decin L, de Koter A, Sana H (2021) Analytic, dust-independent mass-loss rates for red supergiant winds initiated by turbulent pressure. A&A646:A180. https://doi.org/10.1051/0004-6361/202039224. arXiv:2101.03070 [astro-ph.SR]
  * Kervella et al. (2016) Kervella P, Lagadec E, Montargès M, et al. (2016) The close circumstellar environment of Betelgeuse. III. SPHERE/ZIMPOL imaging polarimetry in the visible. A&A 585:A28. https://doi.org/10.1051/0004-6361/201527134. arXiv:1511.04451 [astro-ph.SR]
  * Kervella et al. (2018) Kervella P, Decin L, Richards AMS, et al. (2018) The close circumstellar environment of Betelgeuse. V. Rotation velocity and molecular envelope properties from ALMA. A&A 609:A67. https://doi.org/10.1051/0004-6361/201731761. arXiv:1711.07983 [astro-ph.SR]
  * Kozyreva et al. (2020) Kozyreva A, Nakar E, Waldman R, Blinnikov S, Baklanov P (2020) Shock breakouts from red supergiants: analytical and numerical predictions. MNRAS494(3):3927–3936. https://doi.org/10.1093/mnras/staa924. arXiv:2003.14097 [astro-ph.HE]
  * Kozyreva et al. (2022) Kozyreva A, Klencki J, Filippenko AV, et al. (2022) The Circumstellar Material around the Type IIP SN 2021yja. ApJ934(2):L31. https://doi.org/10.3847/2041-8213/ac835a. arXiv:2207.09976 [astro-ph.HE]
  * Kravchenko et al. (2018) Kravchenko K, Van Eck S, Chiavassa A, et al. (2018) Tomography of cool giant and supergiant star atmospheres. I. Validation of the method. A&A 610:A29. https://doi.org/10.1051/0004-6361/201731530. arXiv:1711.08327 [astro-ph.SR]
  * Kravchenko et al. (2019) Kravchenko K, Chiavassa A, Van Eck S, et al. (2019) Tomography of cool giant and supergiant star atmospheres. II. Signature of convection in the atmosphere of the red supergiant star $\mu$ Cep. A&A 632:A28. https://doi.org/10.1051/0004-6361/201935809. arXiv:1910.04657 [astro-ph.SR]
  * Kravchenko et al. (2021) Kravchenko K, Jorissen A, Van Eck S, et al. (2021) Atmosphere of Betelgeuse before and during the Great Dimming event revealed by tomography. A&A 650:L17. https://doi.org/10.1051/0004-6361/202039801. arXiv:2104.08105 [astro-ph.SR]
  * Langer (2012) Langer N (2012) Presupernova Evolution of Massive Single and Binary Stars. A&AR 50:107–164. https://doi.org/10.1146/annurev-astro-081811-125534. arXiv:1206.5443 [astro-ph.SR]
  * Lèbre et al. (2014) Lèbre A, Aurière M, Fabas N, et al. (2014) Search for surface magnetic fields in Mira stars. First detection in $\chi$ Cygni. A&A 561:A85. https://doi.org/10.1051/0004-6361/201322826. arXiv:1310.4379 [astro-ph.SR]
  * Levesque et al. (2005) Levesque EM, Massey P, Olsen KAG, et al. (2005) The Effective Temperature Scale of Galactic Red Supergiants: Cool, but Not As Cool As We Thought. ApJ628(2):973–985. https://doi.org/10.1086/430901. arXiv:astro-ph/0504337 [astro-ph]
  * Liljegren et al. (2018) Liljegren S, Höfner S, Freytag B, Bladh S (2018) Atmospheres and wind properties of non-spherical AGB stars. A&A 619:A47. https://doi.org/10.1051/0004-6361/201833203. arXiv:1808.05043 [astro-ph.SR]
  * López Ariste et al. (2019) López Ariste A, Tessore B, Carlín ES, et al. (2019) Asymmetric shocks in $\chi$ Cygni observed with linear spectropolarimetry. A&A632:A30. https://doi.org/10.1051/0004-6361/201936189. arXiv:1909.03720 [astro-ph.SR]
  * López Ariste et al. (2022) López Ariste A, Georgiev S, Mathias P, et al. (2022) Three-dimensional imaging of convective cells in the photosphere of Betelgeuse. A&A661:A91. https://doi.org/10.1051/0004-6361/202142271. arXiv:2202.12011 [astro-ph.SR]
  * Ludwig et al. (1999) Ludwig HG, Freytag B, Steffen M (1999) A calibration of the mixing-length for solar-type stars based on hydrodynamical simulations. I. Methodical aspects and results for solar metallicity. A&A346:111–124. https://doi.org/10.48550/arXiv.astro-ph/9811179. arXiv:astro-ph/9811179 [astro-ph]
  * Ludwig et al. (2009) Ludwig HG, Caffau E, Steffen M, et al. (2009) The CIFIST 3D model atmosphere grid. Mem. Soc. Astron. Italiana80:711. https://doi.org/10.48550/arXiv.0908.4496. arXiv:0908.4496 [astro-ph.SR]
  * Ma et al. (2023) Ma JZ, Chiavassa A, de Mink SE, et al. (2023) Is Betelgeuse really rotating? Synthetic ALMA observations of large-scale convection in 3D simulations of Red Supergiants. arXiv e-prints arXiv:2311.16885. https://doi.org/10.48550/arXiv.2311.16885. arXiv:2311.16885 [astro-ph.SR]
  * Magain (1986) Magain P (1986) Contribution functions and the depths of formation of spectral lines. A&A163(1-2):135–139
  * Magg et al. (2022) Magg E, Bergemann M, Serenelli A, et al. (2022) Observational constraints on the origin of the elements. IV. Standard composition of the Sun. A&A661:A140. https://doi.org/10.1051/0004-6361/202142971. arXiv:2203.02255 [astro-ph.SR]
  * Magic and Weiss (2016) Magic Z, Weiss A (2016) Surface-effect corrections for the solar model. A&A592:A24. https://doi.org/10.1051/0004-6361/201527690. arXiv:1606.01030 [astro-ph.SR]
  * Magic et al. (2013) Magic Z, Collet R, Asplund M, et al. (2013) The Stagger-grid: A grid of 3D stellar atmosphere models. I. Methods and general properties. A&A557:A26. https://doi.org/10.1051/0004-6361/201321274. arXiv:1302.2621 [astro-ph.SR]
  * Magic et al. (2015a) Magic Z, Chiavassa A, Collet R, Asplund M (2015a) The Stagger-grid: A grid of 3D stellar atmosphere models. IV. Limb darkening coefficients. A&A573:A90. https://doi.org/10.1051/0004-6361/201423804. arXiv:1403.3487 [astro-ph.SR]
  * Magic et al. (2015b) Magic Z, Weiss A, Asplund M (2015b) The Stagger-grid: A grid of 3D stellar atmosphere models. III. The relation to mixing length convection theory. A&A573:A89. https://doi.org/10.1051/0004-6361/201423760. arXiv:1403.1062 [astro-ph.SR]
  * Massey and Olsen (2003) Massey P, Olsen KAG (2003) The Evolution of Massive Stars. I. Red Supergiants in the Magellanic Clouds. AJ126(6):2867–2886. https://doi.org/10.1086/379558. arXiv:astro-ph/0309272 [astro-ph]
  * Mathias et al. (2018) Mathias P, Aurière M, López Ariste A, et al. (2018) Evolution of the magnetic field of Betelgeuse from 2009-2017. A&A 615:A116. https://doi.org/10.1051/0004-6361/201732542. arXiv:1804.01831 [astro-ph.SR]
  * Merrill (1955) Merrill PW (1955) Red Stars. PASP67(397):199. https://doi.org/10.1086/126804
  * Meynet and Maeder (1997) Meynet G, Maeder A (1997) Stellar evolution with rotation. I. The computational method and the inhibiting effect of the $\mu$-gradient. A&A321:465–476
  * Meynet et al. (2015) Meynet G, Chomienne V, Ekström S, et al. (2015) Impact of mass-loss on the evolution and pre-supernova properties of red supergiants. A&A575:A60. https://doi.org/10.1051/0004-6361/201424671. arXiv:1410.8721 [astro-ph.SR]
  * Montargès et al. (2018) Montargès M, Norris R, Chiavassa A, et al. (2018) The convective photosphere of the red supergiant CE Tauri. I. VLTI/PIONIER H-band interferometric imaging. A&A 614:A12. https://doi.org/10.1051/0004-6361/201731471. arXiv:1802.06086 [astro-ph.SR]
  * Montargès et al. (2021) Montargès M, Cannon E, Lagadec E, et al. (2021) A dusty veil shading Betelgeuse during its Great Dimming. Nature 594(7863):365–368. https://doi.org/10.1038/s41586-021-03546-8
  * Morag et al. (2023) Morag J, Sapir N, Waxman E (2023) Shock cooling emission from explosions of red supergiants - I. A numerically calibrated analytic model. MNRAS522(2):2764–2776. https://doi.org/10.1093/mnras/stad899. arXiv:2207.06179 [astro-ph.HE]
  * Mosumgaard et al. (2018) Mosumgaard JR, Ball WH, Silva Aguirre V, Weiss A, Christensen-Dalsgaard J (2018) Stellar models with calibrated convection and temperature stratification from 3D hydrodynamics simulations. MNRAS478(4):5650–5659. https://doi.org/10.1093/mnras/sty1442. arXiv:1806.00020 [astro-ph.SR]
  * Nakar and Sari (2010) Nakar E, Sari R (2010) Early Supernovae Light Curves Following the Shock Breakout. ApJ725(1):904–921. https://doi.org/10.1088/0004-637X/725/1/904. arXiv:1004.2496 [astro-ph.HE]
  * Nordlund (1982) Nordlund A (1982) Numerical simulations of the solar granulation. I. Basic equations and methods. A&A107:1–10
  * Nordlund et al. (2009) Nordlund Å, Stein RF, Asplund M (2009) Solar Surface Convection. Living Reviews in Solar Physics 6(1):2. https://doi.org/10.12942/lrsp-2009-2
  * Nordlund et al. (2018) Nordlund Å, Ramsey JP, Popovas A, Küffmeier M (2018) DISPATCH: a numerical simulation framework for the exa-scale era - I. Fundamentals. MNRAS477(1):624–638. https://doi.org/10.1093/mnras/sty599. arXiv:1705.10774 [astro-ph.IM]
  * Ohnaka et al. (2017) Ohnaka K, Weigelt G, Hofmann KH (2017) Vigorous atmospheric motion in the red supergiant star Antares. Nature 548(7667):310–312. https://doi.org/10.1038/nature23445. arXiv:1708.06372 [astro-ph.SR]
  * Paladini et al. (2018a) Paladini C, Baron F, Jorissen A, et al. (2018a) Large granulation cells on the surface of the giant star $\pi$1 Gruis. Nature 553(7688):310–312. https://doi.org/10.1038/nature25001
  * Paladini et al. (2018b) Paladini C, Baron F, Jorissen A, et al. (2018b) Constraining Convection in Evolved Stars with the VLTI. The Messenger 172:24–26. https://doi.org/10.18727/0722-6691/5077
  * Paxton et al. (2011) Paxton B, Bildsten L, Dotter A, et al. (2011) Modules for Experiments in Stellar Astrophysics (MESA). ApJS192:3. https://doi.org/10.1088/0067-0049/192/1/3. arXiv:1009.1622 [astro-ph.SR]
  * Paxton et al. (2013) Paxton B, Cantiello M, Arras P, et al. (2013) Modules for Experiments in Stellar Astrophysics (MESA): Planets, Oscillations, Rotation, and Massive Stars. ApJS208:4. https://doi.org/10.1088/0067-0049/208/1/4. arXiv:1301.0319 [astro-ph.SR]
  * Paxton et al. (2015) Paxton B, Marchant P, Schwab J, et al. (2015) Modules for Experiments in Stellar Astrophysics (MESA): Binaries, Pulsations, and Explosions. ApJS220:15. https://doi.org/10.1088/0067-0049/220/1/15. arXiv:1506.03146 [astro-ph.SR]
  * Paxton et al. (2018) Paxton B, Schwab J, Bauer EB, et al. (2018) Modules for Experiments in Stellar Astrophysics (MESA): Convective Boundaries, Element Diffusion, and Massive Star Explosions. ApJS234:34. https://doi.org/10.3847/1538-4365/aaa5a8. arXiv:1710.08424 [astro-ph.SR]
  * Paxton et al. (2019) Paxton B, Smolec R, Schwab J, et al. (2019) Modules for Experiments in Stellar Astrophysics (MESA): Pulsating Variable Stars, Rotation, Convective Boundaries, and Energy Conservation. ApJS243:10. https://doi.org/10.3847/1538-4365/ab2241. arXiv:1903.01426 [astro-ph.SR]
  * Plez (2012) Plez B (2012) Turbospectrum: Code for spectral synthesis. Astrophysics Source Code Library, record ascl:1205.004. ascl:1205.004
  * Rabinak and Waxman (2011) Rabinak I, Waxman E (2011) The Early UV/Optical Emission from Core-collapse Supernovae. ApJ728(1):63. https://doi.org/10.1088/0004-637X/728/1/63. arXiv:1002.3414 [astro-ph.HE]
  * Riva and Steiner (2022) Riva F, Steiner O (2022) Methodology for estimating the magnetic Prandtl number and application to solar surface small-scale dynamo simulations. A&A660:A115. https://doi.org/10.1051/0004-6361/202142644. arXiv:2202.12115 [astro-ph.SR]
  * Roe (1986) Roe PL (1986) Discrete Models for the Numerical Analysis of Time-Dependent Multidimensional Gas Dynamics. Journal of Computational Physics 63(2):458–476. https://doi.org/10.1016/0021-9991(86)90204-4
  * Sagiv et al. (2014) Sagiv I, Gal-Yam A, Ofek EO, et al. (2014) Science with a Wide-field UV Transient Explorer. AJ147(4):79. https://doi.org/10.1088/0004-6256/147/4/79. arXiv:1303.6194 [astro-ph.CO]
  * Salaris and Cassisi (2015) Salaris M, Cassisi S (2015) Stellar models with mixing length and T($\tau$) relations calibrated on 3D convection simulations. A&A577:A60. https://doi.org/10.1051/0004-6361/201525812. arXiv:1503.04582 [astro-ph.SR]
  * Sapir and Waxman (2017) Sapir N, Waxman E (2017) UV/Optical Emission from the Expanding Envelopes of Type II Supernovae. ApJ838(2):130. https://doi.org/10.3847/1538-4357/aa64df. arXiv:1607.03700 [astro-ph.HE]
  * Sapir et al. (2011) Sapir N, Katz B, Waxman E (2011) Non-relativistic Radiation-mediated Shock Breakouts. I. Exact Bolometric Planar Breakout Solutions. ApJ742(1):36. https://doi.org/10.1088/0004-637X/742/1/36. arXiv:1103.5075 [astro-ph.HE]
  * Sapir et al. (2013) Sapir N, Katz B, Waxman E (2013) Non-relativistic Radiation Mediated Shock Breakouts. III. Spectral Properties of Supernova Shock Breakout. ApJ774(1):79. https://doi.org/10.1088/0004-637X/774/1/79. arXiv:1304.6428 [astro-ph.HE]
  * Schawinski et al. (2008) Schawinski K, Justham S, Wolf C, et al. (2008) Supernova Shock Breakout from a Red Supergiant. Science 321(5886):223–226. https://doi.org/10.1126/science.1160456. arXiv:0803.3596 [astro-ph]
  * Schultz et al. (2022) Schultz WC, Bildsten L, Jiang YF (2022) Stochastic Low-frequency Variability in Three-dimensional Radiation Hydrodynamical Models of Massive Star Envelopes. ApJ924(1):L11. https://doi.org/10.3847/2041-8213/ac441f. arXiv:2110.13944 [astro-ph.SR]
  * Schultz et al. (2023) Schultz WC, Bildsten L, Jiang YF (2023) Turbulence-supported Massive Star Envelopes. ApJ951(2):L42. https://doi.org/10.3847/2041-8213/acdf50. arXiv:2306.08034 [astro-ph.SR]
  * Schwarzschild (1952) Schwarzschild M (1952) Shock waves in the atmosphere of pul- sating stars. Transactions of the International Astronomical Union 8:811–812
  * Shussman et al. (2016) Shussman T, Waldman R, Nakar E (2016) Type II supernovae Early Light Curves. arXiv e-prints arXiv:1610.05323 [astro-ph.HE]
  * Shvartzvald et al. (2023) Shvartzvald Y, Waxman E, Gal-Yam A, et al. (2023) ULTRASAT: A wide-field time-domain UV space telescope. arXiv e-prints arXiv:2304.14482. https://doi.org/10.48550/arXiv.2304.14482. arXiv:2304.14482 [astro-ph.IM]
  * Sonoi et al. (2019) Sonoi T, Ludwig HG, Dupret MA, et al. (2019) Calibration of mixing-length parameter $\alpha$ for MLT and FST models by matching with CO5BOLD models. A&A621:A84. https://doi.org/10.1051/0004-6361/201833495. arXiv:1811.05229 [astro-ph.SR]
  * Stein and Nordlund (1998) Stein RF, Nordlund Å (1998) Simulations of Solar Granulation. I. General Properties. ApJ499(2):914–933. https://doi.org/10.1086/305678
  * Stone et al. (2020) Stone JM, Tomida K, White CJ, Felker KG (2020) The Athena++ Adaptive Mesh Refinement Framework: Design and Magnetohydrodynamic Solvers. ApJS249(1):4. https://doi.org/10.3847/1538-4365/ab929b. arXiv:2005.06651 [astro-ph.IM]
  * Stothers and Chin (1995) Stothers RB, Chin CW (1995) Tests of Two Convection Theories for Red Giant and Red Supergiant Envelopes. ApJ440:297. https://doi.org/10.1086/175270
  * Strang (1968) Strang G (1968) On the Construction and Comparison of Difference Schemes. SIAM Journal on Numerical Analysis 5(3):506–517. https://doi.org/10.1137/0705041
  * Tessore et al. (2017) Tessore B, Lèbre A, Morin J, et al. (2017) Measuring surface magnetic fields of red supergiant stars. A&A 603:A129. https://doi.org/10.1051/0004-6361/201730473. arXiv:1704.07761 [astro-ph.SR]
  * Trampedach et al. (2013) Trampedach R, Asplund M, Collet R, Nordlund Å, Stein RF (2013) A Grid of Three-dimensional Stellar Atmosphere Models of Solar Metallicity. I. General Properties, Granulation, and Atmospheric Expansion. ApJ769(1):18. https://doi.org/10.1088/0004-637X/769/1/18. arXiv:1303.1780 [astro-ph.SR]
  * Trampedach et al. (2014) Trampedach R, Stein RF, Christensen-Dalsgaard J, Nordlund Å, Asplund M (2014) Improvements to stellar structure models, based on a grid of 3D convection simulations - II. Calibrating the mixing-length formulation. MNRAS445(4):4366–4384. https://doi.org/10.1093/mnras/stu2084. arXiv:1410.1559 [astro-ph.SR]
  * Wheeler and Chatzopoulos (2023) Wheeler JC, Chatzopoulos E (2023) Betelgeuse: a review. Astronomy and Geophysics 64(3):3.11–3.27. https://doi.org/10.1093/astrogeo/atad020. arXiv:2306.09449 [astro-ph.SR]
  * Wood et al. (1992) Wood PR, Whiteoak JB, Hughes SMG, et al. (1992) OH/IR Stars in the Magellanic Clouds. ApJ 397:552. https://doi.org/10.1086/171812
  * Yoon and Cantiello (2010) Yoon SC, Cantiello M (2010) Evolution of Massive Stars with Pulsation-driven Superwinds During the Red Supergiant Phase. ApJ717(1):L62–L65. https://doi.org/10.1088/2041-8205/717/1/L62. arXiv:1005.4925 [astro-ph.SR]