# Is Betelgeuse really rotating?
Synthetic ALMA observations of large-scale convection in 3D simulations of Red
Supergiants

Jing-Ze Ma (马竟泽) Max Planck Institute for Astrophysics, Karl-Schwarzschild-
Str. 1, 85748 Garching, Germany<EMAIL_ADDRESS>Andrea Chiavassa
Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, Lagrange, CS
34229, Nice, France  Max Planck Institute for Astrophysics, Karl-
Schwarzschild-Str. 1, 85748 Garching, Germany  Selma E. de Mink Max Planck
Institute for Astrophysics, Karl-Schwarzschild-Str. 1, 85748 Garching, Germany
Anton Pannekoek Institute for Astronomy and GRAPPA, University of Amsterdam,
NL-1090 GE Amsterdam, The Netherlands  Ruggero Valli Max Planck Institute for
Astrophysics, Karl-Schwarzschild-Str. 1, 85748 Garching, Germany  Stephen
Justham Max Planck Institute for Astrophysics, Karl-Schwarzschild-Str. 1,
85748 Garching, Germany  School of Astronomy & Space Science, University of
the Chinese Academy of Sciences, Beĳing 100012, China  National Astronomical
Observatories, Chinese Academy of Sciences, Beĳing 100012, China  Bernd
Freytag Theoretical Astrophysics, Department of Physics and Astronomy, Uppsala
University, Box 516, SE-751 20 Uppsala, Sweden

###### Abstract

The evolved stages of massive stars are poorly understood, but invaluable
constraints can be derived from spatially resolved observations of nearby red
supergiants, such as Betelgeuse. ALMA observations of Betelgeuse showing a
dipolar velocity field have been interpreted as evidence for a rotation rate
of $v\sin i\sim 5\,\mathrm{km\,s^{-1}}$. This is two orders of magnitude
larger than predicted by single-star evolution, leading to the suggestion that
Betelgeuse is a binary merger product. We propose instead that the velocity
field could be due to large-scale convective motions. The resulting surface
velocity maps can sometimes be mistaken for rotation, especially when the
turbulent motions are only partially resolved, as is the case for the current
ALMA beam. We support this claim with 3D CO5BOLD simulations of _non-rotating_
red supergiants post-processed to predict synthetic ALMA images and SiO
spectra to compare with observed radial velocity maps. Our simulations show a
$\sim 50\%$ chance to be interpreted as evidence for a rotation rate as high
as claimed for Betelgeuse. We conclude that we need at least another ALMA
observation to firmly establish whether Betelgeuse is indeed rapidly rotating.
Such observations would also provide insight into the role of angular momentum
and binary interaction in the late evolutionary stages. The data will further
probe the structure and complex physical processes in the atmospheres of red
supergiants, which are immediate progenitors of supernovae and are believed to
be essential in the formation of gravitational wave sources.

††software:  CO5BOLD (Freytag et al., 2012), FASTCHEM2 (Stock et al., 2018,
2022), MARCS (Gustafsson et al., 2008), MAGRITTE (De Ceuster et al., 2020a, b,
2022), MESA (Paxton et al., 2011, 2013, 2015, 2018, 2019; Jermyn et al.,
2023), Astropy (Astropy Collaboration et al., 2013, 2018, 2022), NumPy (Harris
et al., 2020), SciPy (Virtanen et al., 2020), Matplotlib (Hunter, 2007),
Jupyter (Kluyver et al., 2016)

##  1 Introduction

Cool evolved stars are not expected to be rotating fast, at least not at their
surfaces. As the stars evolve, their envelopes expand by one to two orders of
magnitude. The outer layers thus slow down as a result of angular momentum
conservation and may be further reduced by, e.g. mass loss due to stellar
winds (e.g. Maeder & Meynet, 2000; Smith, 2014), possibly inward convective
transport of angular momentum (e.g. Brun & Toomre, 2002; Brun & Palacios,
2009), and magnetic braking (Mestel, 1968). The theory of single star
evolution therefore predicts slow surface rotation rates, less than about
$1\,\mathrm{km\,s^{-1}}$ for stars at the tip of the red giant branch (e.g.
Privitera et al., 2016a), and less than about $0.1\,\mathrm{km\,s^{-1}}$ for
red supergiants (RSGs; Wheeler et al., 2017; Chatzopoulos et al., 2020), which
are the cool giant descendants of massive stars.

Despite theoretical expectations, cool stars with rotation rates exceeding
these predictions have been observed across the Hertzsprung-Russell diagram.
These include several hundred red giants, about $1\%$ of the total population
of red giants (e.g. Patton et al., 2023, and references therein), and a few
asymptotic giant branch (AGB) stars (Barnbaum et al., 1995; Vlemmings et al.,
2018; Brunner et al., 2019; Nhung et al., 2021, 2023). For RSGs, so far only
one has been claimed to rotate rapidly: $\alpha$ Orionis, better known as
Betelgeuse (Uitenbroek et al., 1998; Harper & Brown, 2006; Kervella et al.,
2018), which recently has drawn wide attention after the sudden Great Dimming
(Guinan et al., 2019) and subsequent re-brightening (Guinan et al., 2020).

Betelgeuse, being one of the closest RSGs to Earth, is one of the few stars
that can be spatially resolved and has therefore been a target of
interferometric studies for over a century (Michelson & Pease, 1921).
Recently, the Atacama Large Millimeter/submillimeter Array (ALMA) provided
unprecedented maps of the molecular envelope (Kervella et al., 2018, hereafter
K18; right-hand panels of our Figure 2). The surface radial velocity map shows
a remarkably clear dipolar structure: half of the visible hemisphere of the
star shows a blue shift and the other half shows a red shift of several
$\mathrm{km\,s^{-1}}$.

A natural explanation of such a dipolar velocity field is stellar rotation, as
noted by K18. They inferred a projected equatorial velocity of $v\sin
i=5.47\pm 0.25\,\mathrm{km\,s^{-1}}$. They compared the results with earlier
measurements using the Hubble Space Telescope (HST) probing the chromosphere
(Uitenbroek et al., 1998; Harper & Brown, 2006) and argued that both the ALMA
and HST data are consistent with the interpretation that Betelgeuse is fast-
rotating. The fast rotation rate inferred for Betelgeuse is surprising in
light of the predictions of single-star models, as illustrated in Figure 1 and
Appendix A for details.

Binary star evolution has been proposed as an explanation for Betelgeuse’s
high rotation rate, in particular the merger with a lower mass companion
(Wheeler et al., 2017; Chatzopoulos et al., 2020; Sullivan et al., 2020;
Shiber et al., 2023). This may seem like an exotic explanation, but massive
stars often interact with close companions (Sana et al., 2012). As a
consequence, stellar mergers are expected to be common (de Mink et al., 2014;
Kochanek et al., 2014). Zapartas et al. (2019) estimated that as many as one-
third of RSGs experience a stellar merger before they reach core collapse. Rui
& Fuller (2021) identified two dozen red giants that are possible merger
products, based on their asteroseismological signatures. For red giants, the
engulfment of planets has also been proposed as an explanation for their rapid
rotation (Carlberg et al., 2012; Privitera et al., 2016b; Gaulme et al., 2020;
Lau et al., 2022).

Figure 1: The equatorial rotational velocity inferred for Betelgeuse ( K18, )
is 2–3 orders of magnitude higher than single-star models predict. To
illustrate this point, we show MESA predictions for the rotational velocities
of massive stars ($\sim 10-20\,M_{\odot}$) with typical initial rotation rates
as a function of their fractional age, i.e. time divided by their main
sequence lifetime. The steep drop in rotation rates marks the rapid expansion
as stars transition from the core H-burning phase (blue background) to the
core He-burning RSGs (red background). See Appendix A for details.

Establishing whether Betelgeuse is indeed rotating, is of vital importance to
better understand its evolutionary history, the possible role of binary
interaction, and the physics of the evolved stages of massive stars in general
(see Wheeler & Chatzopoulos, 2023, for a review). Unfortunately, accurate and
reliable measurements of rotation rates for red (super)giants are challenging.

The first complicating factor concerns the high velocities expected for
convective flows at the photosphere. These may be as high as
$20\,\mathrm{km\,s^{-1}}$, as shown in different 3D simulations (Kravchenko et
al., 2019; Antoni & Quataert, 2022; Goldberg et al., 2022) as well as
spectroscopic (Lobel & Dupree, 2000; Josselin & Plez, 2007) and optical
interferometric observations (Ohnaka et al., 2009, 2011, 2013, 2017). This is
two orders of magnitude larger than the predicted rotational velocities and
four times larger than the rotational velocity inferred for Betelgeuse by K18.
How the turbulent velocity field affects the measurement of the rotation rate
is not yet well understood.

A second complication concerns the large sizes expected for the convective
cells at the surface, which may span a significant fraction of the radius
(Schwarzschild, 1975). Only a few of them will be present at the surface at
any given time, as also suggested by spectropolarimetric (Ariste et al., 2018,
2022) and optical interferometric observations (e.g. Haubois et al., 2009;
Norris et al., 2021). If, by chance, one very large cell or a group of cells
move toward the observer while others move away, this can result in a dipolar
velocity field even for a _non-rotating_ star.

The central question motivating our current study is: “Can a _non-rotating_
red (super)giant be mistaken to be a rapid rotator?” For our study, we use
existing 3D radiation hydrodynamic simulations of RSGs with properties similar
to Betelgeuse. We develop a new post-processing package to solve the radiative
transfer equations and make direct predictions for ALMA observables that we
compare with observations of Betelgeuse. We quantify how fast a _non-rotating_
star can appear to be rotating and how likely it is to obtain spurious
measurements of high rotation. We conclude that, to firmly establish whether
Betelgeuse is rotating rapidly, additional epochs of ALMA observations are
needed, preferably with higher spatial resolution.

##  2 Method

###  2.1 3D simulations of RSG envelope with CO5BOLD

To assess whether a _non-rotating_ RSG can show a dipolar radial velocity map
in the ALMA band, we need global 3D RSG models that simulate the multi-scale
convection of the full star. So far this is only possible with the CO5BOLD
models (Freytag et al., 2012), since other 3D RSG models do not simulate the
whole $4\pi$ sphere (Goldberg et al., 2022). We use existing 3D radiation-
hydrodynamic simulations of RSG convective envelopes from Ahmad et al. (2023),
computed with the CO5BOLD code. The CO5BOLD RSG simulations have been
extensively used to interpret spectro-photometric, interferometric, and
astrometric observations, especially in the context of Betelgeuse (e.g.
Chiavassa et al., 2009, 2010; Montargès et al., 2014, 2016; Kravchenko et al.,
2021). The code numerically integrates the nonlinear compressible hydrodynamic
equations, coupled with a short-characteristics scheme for radiation transport
(Freytag et al., 2012). The global simulations that we use in this work adopt
the ‘star-in-a-box’ setup by simulating the outer part of the convective
envelope with mass $M_{\mathrm{env}}$ on an equidistant Cartesian grid. The
interior is replaced by an artificial central region providing a luminosity
source with a drag force to damp the velocity. A gravitational potential is
imposed, set by the central mass $M_{\mathrm{pot}}$. For the gravity
experienced by the simulation, the stellar mass is $M_{\mathrm{pot}}$ because
the self-gravity of the envelope is neglected. However, for actual stars, the
total stellar mass would be $M_{\mathrm{pot}}+M_{\mathrm{env}}$. For detailed
descriptions of the general setup, see Chiavassa et al. (2011a), Freytag et
al. (2012) and Chiavassa et al. (2024, submitted).

We use the simulation of a _non-rotating_ star with large surface convective
cells, given by the surface gravity $\log_{10}\,g=-0.419$ in cgs units and
effective temperature $T_{\mathrm{eff}}=3366$ K, close to the observed values
for Betelgeuse (see Table 2 in Appendix D). An alternative simulation, with
smaller surface convective cell sizes but which fits better with the mass and
luminosity inferred for Betelgeuse, is presented and discussed in Appendix D.

###  2.2 Synthetic ALMA images

We post-process the 3D simulations to create synthetic observations for ALMA.
Here we describe the main assumptions. Details and tests for each step are
provided in Appendix C. The post-processing package, animations and data
behind the figures are publicly available at Zenodo: doi:
10.5281/zenodo.10199936.

To directly compare with the SiO line spectra observed by ALMA, we need to
calculate the intensity from SiO emission and absorption. For the chemical
abundances of the relevant species (SiO molecules, the electrons and atomic
H), we use the equilibrium-chemistry code
FASTCHEM2111https://github.com/exoclime/FastChem (Stock et al., 2018, 2022).
FASTCHEM2 has been widely adopted and tested against observations of
exoplanetary atmosphere, in particular for ultra-hot jupiters where the day-
side conditions are similar to cool stars (e.g. Kitzmann et al., 2018, 2023).
For AGB atmospheres which share similar convective properties with RSGs,
chemical equilibrium was shown to be a reasonable approximation (Agúndez et
al., 2020). Here, we present the first application of FASTCHEM2 to stellar
atmospheres.

To synthesize the intensity map, we numerically integrate the radiative
transfer equation on a Cartesian grid. We take into account the continuum
free-free opacity and Doppler-shifted $\mathrm{{}^{28}SiO}$ lines (vibrational
level $\nu=2$, rotational transition $J=8-7$ as observed by ALMA). Detailed
equations and opacity sources can be found in Appendix C.2.

We scale the simulation such that the radio photospheric radius approximately
matches ALMA observations in K18 (For the associated limitations, see Appendix
C.4). We interpolate the intensities onto a $101^{2}$ pixel grid, such that
both the pixel number and the angular coverage are identical to ALMA
observations. We then convolve the intensity map with the ALMA beam (FWHM of
18 mas). We use the convolved continuum intensity map to calculate the radius
of the radio photosphere $R_{\mathrm{star}}=\sqrt{\int
I(r,\theta)r\mathrm{d}r\mathrm{d}\theta/(\pi I_{\mathrm{star}})}$ where the
integral is performed over the 2D image (Section 3.3.1 in K18, ). Here, $I$ is
the continuum intensity as a function of the polar coordinate $(r,\theta)$ on
the 2D image, measured in the spectral window centered at $343.38$ GHz as done
in ALMA observations. Its value at the photocenter is denoted as
$I_{\mathrm{star}}$.

Figure 2:  A direct comparison of a _non-rotating_ RSG simulation with ALMA
observations of Betelgeuse. For this figure, we have chosen an illustrative
example of a simulation snapshot in time and orientation, where the _non-
rotating_ simulation can be easily mistaken as rotation. The left panels show
the unconvolved simulation image (a & d). The middle panels show the convolved
synthetic image (b & e). The right panels show the actual ALMA observations (c
& f) from K18. The ALMA beam size is indicated by the gray circles. The top
row demonstrates the continuum intensity maps. The measured radii
$R_{\mathrm{star}}$ of the radio photosphere are indicated in white circles.
The bottom row shows the radial velocity maps measured from the Doppler shift
of SiO lines. The green arrows indicate the rotational axes from the fit.
Images from the simulation are rotated such that the rotational axis aligns
with the observed axis. Credit for the right two panels: Kervella et al., A&A,
609, A67, 2018, reproduced with permission © ESO. This figure is available as
an animation in the HTML version of the paper as well as on Zenodo: doi:
10.5281/zenodo.10199936 for our 5-year-long mock observations.

###  2.3 Analysis to the radial velocity map

To compare our simulations with the observed radial velocity maps, we follow
the analysis in K18. We use the continuum-subtracted intensities to identify
the line in each pixel, fit a Gaussian profile to the line using least squares
fitting, and obtain the radial velocity shift from the central value of the
Gaussian. The apparent systematic velocity $v_{\mathrm{sys}}$ is obtained from
the radial velocity shift of the integrated line (over the region up to $30$
mas from the center for the absorption line and $30-50$ mas for the emission
line as in Section 3.3.2 in K18).

We measure the apparent $v\sin i$ by fitting the projected radial velocity map
of a rigidly-rotating sphere to the $v_{\mathrm{sys}}$-subtracted radial
velocity map within $R_{\mathrm{star}}$ (Section 3.5 in K18, however, see a
more detailed discussion in Appendix C.3 for the choice of the radius
adopted). Since both the observing time and the orientation of the star with
respect to the observer are arbitrary, we repeat this procedure for six faces
of the Cartesian box and every snapshot of the simulations. Throughout the
main text, we use the radial velocity maps derived from emission lines in
accordance with the actual observations of K18.

##  3 Results

###  3.1 Selected mock observations compared to actual ALMA observations

In Figure 2, we show a selected example of a synthetic ALMA image from our
_non-rotating_ RSG simulation, and compare it with the ALMA observations of
Betelgeuse in K18. The left-hand panels (a & d) show the original simulation.
The middle panels (b & e) show our synthetic observations after convolution
with the ALMA beam. The right-hand panels (c & f) show the actual ALMA
observations ( K18, ). In the top row, we show the simulated and observed
continuum intensity maps. In the bottom row, we show the radial velocity maps
measured from the Doppler shifts of the lines.

The unconvolved intensity map at these wavelengths probes the atmospheric
layers of the star which are highly asymmetric (Figure 1a). The convective
motions of $\pm\,30\,\mathrm{km\,s^{-1}}$ can be seen in the radial velocity
map of the original simulations (Figure 1d). Both the cell size and the peak
convective velocity are consistent with analytical estimates (Appendix B) and
other 3D simulations of RSGs (Kravchenko et al., 2019; Antoni & Quataert,
2022; Goldberg et al., 2022).

The ALMA beam used in the settings by K18 is about $30\%$ of the diameter of
the radio photosphere. This means that any sharp features on the surface will
be smoothed out. The effect of this can be seen when comparing the left and
middle panels in Figure 1. For the intensity map, a single hot spot emerges
(Figure 1b), similar to what is observed by ALMA (Figure 1c). Such a feature
has also been observed in other images of Betelgeuse in the UV (Gilliland &
Dupree, 1996), optical (Young et al., 2000; Haubois et al., 2009; Montargès et
al., 2016), and radio (O’Gorman et al., 2017). For the radial velocity map,
positive and negative radial velocity shifts partially cancel each other after
convolution with the beam size. This blurs the original radial velocity map
with convective motions up to $30\,\mathrm{km\,s^{-1}}$ (Figure 1d) into a map
with radial velocity variations up to $8\,\mathrm{km\,s^{-1}}$ (Figure 1e)
consistent with the ALMA observations. We conclude that the synthetic images
of the intensity map and radial velocity qualitatively match the main features
of the observations. In particular, both synthetic and observed radial
velocity maps (Figure 1e and f) show a dipolar velocity field, with one
hemisphere approaching the observer and the other receding.

K18 interpreted the dipolar velocity field as a sign of rotation. Fitting for
this, assuming rigid rotation, they inferred a $v\sin i$ of $5.47\pm
0.10\,\mathrm{km\,s^{-1}}$, and a residual velocity dispersion of
$1.44\,\mathrm{km\,s^{-1}}$. Following the same procedure for our synthetic
image, we obtain very similar values, see Table 1. The resulting
$\chi^{2}_{\mathrm{red}}$ is larger than observed, but comparable within an
order of magnitude. As shown in Figure 10 in Appendix D, if we have
underestimated the overall smearing effect of the observational pipeline, the
$\chi^{2}_{\mathrm{red}}$ may be vastly decreased.

Table 1: Comparison of fitting parameters of a _non-rotating_ RSG simulation with ALMA observations of Betelgeuse. The radial velocity maps being fitted are shown in Figure 2. A rigidly-rotating model is used to fit the radial velocity maps to obtain the perceived projected rotation rate $v\sin i$ and fitting parameters (standard deviation of the residual $\sigma_{\mathrm{res}}$ and reduced $\chi_{\mathrm{red}}^{2}$). All the parameters are obtained using the method described in Section 2.3.  Inferred parameters | $v\sin i$ | $\sigma_{\mathrm{res}}$ | $\chi^{2}_{\mathrm{red}}$
---|---|---|---
$\mathrm{[km\,s^{-1}]}$ | $\mathrm{[km\,s^{-1}]}$ |
Mock observation | $5.51\pm 0.21$ | $1.68$ | $94.6$
(this work, non-rotating)
ALMA observation | $5.47\pm 0.10$ | $1.44$ | $20.7$
(K18)

Note, however, that our simulation is based on a _non-rotating_ star. The
dipolar velocity field shown in the synthetic image is not the result of
rotation, but the alignment of the group motions of the surface convective
cells. The example shown here illustrates how turbulent motions in the
atmosphere of a _non-rotating_ RSG can, at certain times, mimic the effect of
rotation. This raises the question of how confident we are that Betelgeuse is
rotating and to what extent the expected large-scale convection affects the
rotation measurement.

###  3.2 Probability of inferring rotation from the radial velocity map of a
non-rotating RSG

Figure 3:  Probability distribution of apparent $v\sin i$ resulting from
turbulent motions above the surface of a _non-rotating_ stellar model. The
grey vertical lines indicate the observed values based on SiO and CO lines
(K18), with their uncertainties indicated by the widths and the SiO line we
have modeled shown darker. The blue curve illustrates the kernel-smoothed
probability density function of the inferred $v\sin i$. The $50\%$ probability
interval around the peak is shaded in grey. The orange curve shows the
cumulative distribution function.

In the previous section, we showed an example of how a _non-rotating_
simulation can appear to be rotating as a result of under-resolved convective
motions. We had chosen a particular snapshot in time that illustrates this
point well. In this section, we discuss how often such situations occur. We
present the distribution of inferred $v\sin i$ in our mock observations for
the simulated RSG. To obtain the probability distribution of inferred $v\sin
i$ in a single-epoch observation, we compile 480 inferred $v\sin i$ calculated
for 6 faces and 80 snapshots (uniformly taken across the 5-year time span of
the relaxed simulation) and plot the distribution in Figure 3.

Figure 3 shows that the apparent $v\sin i$ distribution inferred for our _non-
rotating_ models peaks at $\sim 5\,\mathrm{km\,s^{-1}}$. The observed values
for Betelgeuse (grey vertical lines) fall into the $50\%$ probability interval
around the peak. These simulations thus show that it is rather common that
turbulent motions give rise to apparent $v\sin i$ similar to what is observed
for Betelgeuse. We further show that for a _non-rotating_ RSG to be
interpreted as rotating faster than $\sim 5\,\mathrm{km\,s^{-1}}$, the
probability still remains as high as $\sim 40\%$. These numbers depend on the
simulation adopted, and the parameters assumed for Betelgeuse, which are not
well-constrained (Dolan et al., 2016; Joyce et al., 2020). Nevertheless, from
the distribution we expect that $\sim 90\%$ single-epoch observations would
show signs of rotation faster than $\sim 2\,\mathrm{km\,s^{-1}}$. By visual
inspection, the majority of those with inferred $v\sin i\lesssim
2\,\mathrm{km\,s^{-1}}$ do not show an obvious dipolar structure that could be
mistaken for rotation.

Figure 4:  Three interpretations of Betelgeuse’s dipolar radial velocity map
observed by ALMA. The physical pictures and their observational consequences
are explained in Section 4.

##  4 Discussion

In Section 3, we argued that convective motions may be responsible for the
dipolar velocity field observed for Betelgeuse, as an alternative to the
explanation of rapid rotation. The quantitative results we presented are
subject to uncertainties, due to the choices for the model parameters and
limitations of simulations. We discuss this further in Appendix C.4 & D.
Despite the uncertainties, we can use these results to formulate conceptual
physical pictures that can be tested with future observations. In Figure 4 we
illustrate three scenarios that, in principle, can explain the dipolar
velocity map observed for the molecular shell around Betelgeuse.

  1. I.

Rapid rotation: This hypothesis, proposed by K18, states that Betelgeuse is
rapidly rotating and drags the surrounding molecular shell along.

  2. II.

Large convective cells and stochastic effects: In this hypothesis, only very
few convective cells cover the surface leading to stochastic effects (e.g.
Schwarzschild, 1975). At certain times, only two large cells may dominate the
dynamics of the hemisphere. If one cell moves outward and one inward, and the
motions are transported to the molecular shell via waves or shocks, this would
result in a dipolar velocity field.

  3. III.

Smaller convective cells moving coherently: A dipolar velocity field may also
arise when the convective cells are smaller but move semi-coherently in
groups. This physical picture most closely describes what we see in our
simulations, where coherent motion is the result of deeper convective motions
that operate over a length scale comparable to the stellar radius. Other
mechanisms, e.g. non-radial oscillations (Lobel & Dupree, 2001), may also be
able to drive such coherent motions.

Note that, in principle, the rotational velocity field and the turbulent
motions can co-exist in the molecular layer. The three hypotheses we list here
emphasize which of the components _dominate_ the velocity field.

###  4.1 Current observational constraints

All three hypotheses we consider can reproduce the dipolar radial velocity
map. Here we mention further observational constraints that argue in favor or
against the scenarios depicted in Fig. 4.

Support for the rotation hypothesis (I) has been claimed based on HST data
probing the chromosphere of Betelgeuse (Uitenbroek et al., 1998; Harper &
Brown, 2006). They found an upward trend in the radial velocity while scanning
over the surface from northwest to southeast, and interpreted this as being
the result of rotation. The magnitude and direction of the inferred rotation
are consistent with the data taken by K18 almost 20 years later.

This interpretation of the HST data is not without controversy. The radial
velocities vary between observations taken a few months apart (see Fig. 6 in
Harper & Brown, 2006). Other HST data analyzed by Lobel & Dupree (2001)
instead shows evidence for a reversal of velocities and Jadlovský et al.
(2023) found no sign of rotation. These findings are more consistent with
hypotheses (II) and (III), although it is unclear what the effect is of the
great dimming event in the last study.

Evidence for the presence of (large) convective cells can be found in the
asymmetries in the surface intensity map observed by ALMA, but also in earlier
data probing the UV, optical, and radio, as discussed in Sect. 3.
Specifically, the infrared images taken with VLTI/MATISSE, which probe the
bottom part of the molecular layer with a spatial resolution of 4 mas (Drevon
et al., 2024), show multiple resolved structures. Other nearby RSGs such as
Antares, V602 Car and AZ Cyg show similar evidence (Ohnaka et al., 2017;
Climent et al., 2020; Norris et al., 2021).

A further indication for turbulent motions at the molecular shell is the
integrated line broadening observed by ALMA. The observed broadening of FWHM
$24\,\mathrm{km\,s^{-1}}$ (Section 3.3.2 in K18, ) is significantly larger
than what can be explained by the claimed rotation rate of
$5.5\,\mathrm{km\,s^{-1}}$. This line width is consistent with what we predict
in our simulations, about $18\,\mathrm{km\,s^{-1}}$ on average. There is thus
evidence for additional broadening, likely originating from convective
motions, consistent with hypotheses (II) and (III).

###  4.2 Future observations to test the hypotheses

The hypotheses we propose have clear predictions that can be tested with
future observations. In particular, additional epochs and higher spatial
resolution are desired to resolve the variable radial velocity map.

Additional epochs will allow us to probe changes in the velocity field over
time. In the case of rapid rotation (hypothesis I), we expect that the radial
velocity map will not change significantly. Both the magnitude of the rotation
and alignment of the spin axis should be close to values obtained in 2015 by
K18. Instead, if the velocity field is dominated by turbulent motions, we
expect the surface to readjust. Convective cells at the surface (hypothesis
II) are expected to change on timescales of a few months (see e.g. Montargès
et al., 2018; Norris et al., 2021, for RSG CE Tauri and AZ Cyg). Convective
motions in the deep layers (hypothesis III) may take years, see Appendix B for
an order of magnitude estimate. We expect the surface velocity field to
significantly change on these timescales. The field may still be dipolar, but
likely with a different orientation, or not display a dipolar feature.

Based on the 3D simulations we presented in this work, we estimate that in a
single-epoch observation, there is a $50\%$ chance that a _non-rotating_ RSG
appears to rotate as fast as observed (Figure 3). Assuming that there is no
preferred direction and allowing for 20-degree uncertainty in the rotational
axis, we expect a $20/360\times 50\%\simeq 2\%$ chance that a second epoch for
a non-rotating RSG will give a similar radial velocity map as observed in 2015
by coincidence. Therefore, this idealized example indicates that one
additional single-epoch observation would be sufficient to test the rotation
hypothesis at a confidence level equivalent to two sigma.

Observations of the radial velocity map at increased spatial resolution will
also help distinguish the hypotheses we have outlined here. In Figure 10 in
the Appendix D, we present mock observations with different spatial
resolutions. We expect that an increase of a factor of two in spatial
resolution will already marginally resolve the convective structure both in
the radial velocity map and the intensity map, and be able to distinguish
between all three scenarios.

An increase by a factor of two should be feasible. For ALMA, higher resolution
can in principle already be achieved by going to higher frequencies. Asaki et
al. (2023) reported to have achieved a spatial resolution of 5 mas for a
carbon-rich AGB star, which is three times better. When ALMA upgrades its
current 16 km baseline to 32 km, we can expect an (additional) improvement of
a factor of two (Carpenter et al., 2020).

##  5 Conclusions

In this work, we investigate whether large-scale convection can affect the
rotation measurement of Betelgeuse. We generate synthetic ALMA images from 3D
CO5BOLD simulations of _non-rotating_ RSGs, and compare them to the actual
ALMA observations of Betelgeuse (K18). Our conclusions are summarised as
follows:

  * $\bullet$

Large-scale motions in the RSG atmospheres generated by convection can be
mistaken for rotation in interferometric observations. Both the convective
cell size of the RSGs and the beam size used in ALMA observations are a large
fraction of the stellar radius. Therefore, large-scale convection can be
blurred as a dipolar feature in the radial velocity map that resembles
rotation. This may apply to other cool stars, e.g. the $\sim 1\ M_{\odot}$ AGB
star R Dor, which is reported to be rotating at $v\sin i\sim
1\,\mathrm{km\,s^{-1}}$ with ALMA (Figure 1 in Vlemmings et al., 2018) but
only displays turbulent motions in higher-resolution VLTI/AMBER observations
probing similar heights above the photosphere (Figure 7 in Ohnaka et al.,
2019).

  * $\bullet$

Future interferometric observations, e.g. with ALMA, are needed to assess the
rotation of Betelgeuse. Our simulations suggest that another single-epoch
observation of Betelgeuse with ALMA is sufficient to confirm if the observed
maps show signs of rotation. Multi-epoch observations and higher spatial
resolutions are desired for further constraints.

  * $\bullet$

The post-processing package developed in this work can be applied to other
forward modelings from 3D simulations of (binary) cool stars to synthetic
radio spectra (e.g. Ramstedt et al., 2017; Doan et al., 2020; De Ceuster et
al., 2022). This is particularly timely with the ongoing ALMA programs such as
DEATHSTAR (Ramstedt et al., 2020; Andriantsaralaza et al., 2021) and ATOMIUM
(Decin et al., 2020; Gottlieb et al., 2022; Montargès et al., 2023; Decin et
al., 2023) that advance our understanding of chemistry, dust formation,
planetary nebulae, binary interaction, and mass loss of these cool stars.

  * $\bullet$

Regardless of whether Betelgeuse is rapidly rotating, more efforts are needed
in both theoretical and observational aspects of RSGs. If Betelgeuse is indeed
rapidly rotating, understanding the consequences of convection may well still
be important for accurately interpreting the observational signatures of
rotation. The stellar merger scenario for the rapid rotation (Wheeler et al.,
2017; Chatzopoulos et al., 2020; Sullivan et al., 2020; Shiber et al., 2023)
would demand further studies in the context of Betelgeuse’s other properties,
e.g. its runaway nature (Harper et al., 2008; Decin et al., 2012). If, on the
other hand, future observations show evidence of turbulent motions, it may be
possible for ALMA and other interferometers to trace the velocity field in
different heights across the RSG atmosphere. Such observations would enable us
to investigate the possible connections between pulsation, convection, and
wind-launching mechanisms (Yoon & Cantiello, 2010; Kee et al., 2021).

We acknowledge helpful discussions with Pierre Kervella, Norbert Langer and
Falk Herwig. A.C. acknowledges support from the French National Research
Agency (ANR) funded project PEPPER (ANR-20-CE31-0002). B.F. has received
funding from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation programme (Grant agreement No. 883867,
project EXWINGS). The CO5BOLD simulations were enabled by resources provided
by the Swedish National Infrastructure for Computing (SNIC) and PSMN in Lyon,
France. This research project was partly conducted using computational
resources (and/or scientific computing services) at the Max-Planck Computing &
Data Facility. This research was supported by the Munich Institute for Astro-,
Particle and BioPhysics (MIAPbP) which is funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany´s
Excellence Strategy – EXC-2094 – 390783311.

##  Appendix A Single star 1D models

The fast rotation rate inferred for Betelgeuse by K18 is surprising in light
of the predictions of single-star models (Wheeler et al., 2017; Chatzopoulos
et al., 2020). We illustrated this in Figure 1 mentioned briefly in the main
text. Here we provide additional details.

We show predictions for massive star evolutionary models. During their main-
sequence evolution (core H burning, blue shading), massive stars rotate up to
a few hundred $\mathrm{km\,s^{-1}}$, but once they have expanded to become
RSGs (core He burning, red shading) their rotation rates drop to 0.1-0.001
$\mathrm{km\,s^{-1}}$. In general, the RSGs with higher initial masses rotate
slower due to higher angular momentum loss via stronger stellar wind.

The inferred rotation rate for Betelgeuse is 2-3 orders of magnitude larger
than what is predicted, violating the expectations for single-star evolution.
The figure shows the measurement for Betelgeuse by K18 as a red box. The
bottom of the box corresponds to an assumed orientation where we see the star
equator on ($i=90^{\circ}$). The center of the box assumes $\langle\sin
i\rangle=\pi/4$, which corresponds to the average for random orientations. The
top of the box assumes $i=20^{\circ}$, which is the value originally suggested
by Uitenbroek et al. (1998) based on the location of the hot spot.

The models shown are computed with 1D stellar evolution code MESA v22.05.1
(Paxton et al., 2011, 2013, 2015, 2018, 2019; Jermyn et al., 2023), with
initial masses of $12,\,15,\,18,\,21\,M_{\odot}$ and initial rotational
velocities of $100,\,200,\,300\,\mathrm{km\,s^{-1}}$ at metallicity $Z=0.02$.
Convection is modelled according to mixing-length theory with $\alpha_{\rm
MLT}=1.5$, plus an exponential overshoot scheme (Herwig, 2000) with
coefficients $f_{\rm{ov}}=0.0415$ and $f_{0,\rm{ov}}=0.008$ (Brott et al.,
2011). The convective boundary is determined with the Ledoux (1947) criterion,
and we assume efficient semiconvection (Langer et al., 1983, 1985). During the
main sequence, wind mass loss is accounted for as in Vink et al. (2001), while
we switch to de Jager et al. (1988) when the effective temperature drops below
$10^{4}\,\rm{K}$. We include the following mechanisms of rotational mixing and
angular momentum transport (see Heger et al., 2000; Paxton et al., 2013):
dynamical shear instability, secular shear instability, Eddington-Sweet
circulation, Goldreich-Schubert-Fricke instability and Spruit-Tayler dynamo
(Spruit, 2002; Heger et al., 2005; Paxton et al., 2013).

##  Appendix B Order-of-magnitude estimates

In this section, we estimate the properties of convection of RSGs with order-
of-magnitude analytical arguments. The surface convective cell size is
determined by the effective temperature $T_{\mathrm{eff}}$ and surface gravity
$g$. The convective velocity, on the other hand, is determined by
$T_{\mathrm{eff}}$ alone. These parameters can be estimated on the order-of-
magnitude level, given our understanding of RSGs and measured properties of
Betelgeuse.

Near the photosphere, the horizontal size $d$ of the convective cell can be
estimated from the pressure scale height $H_{p}$ multiplied by a factor
$\alpha$ of order $10$ for RSGs (Tremblay et al., 2013; Paladini et al.,
2018). Therefore, the surface convective cell size is

$d\sim\alpha\frac{k_{\mathrm{B}}T_{\mathrm{eff}}}{\mu m_{\mathrm{u}}g}\sim
66\,\mathrm{R_{\odot}}\left(\frac{\alpha}{10}\right)\left(\frac{T_{\mathrm{eff}}}{3600\,\mathrm{K}}\right)\left(\frac{g}{\mathrm{g\,cm^{-3}}}\right)^{-1}\,,$
(B1)

where $k_{\mathrm{B}}$ is the Boltzmann constant, $m_{\mathrm{u}}$ is the
proton mass, $\mu$ is the mean molecular weight taken to be $0.643$ for near
solar composition. A similar method was also adopted to estimate the size of
granules in 3D CO5BOLD RSG simulations (Chiavassa et al., 2011b). This
estimate indicates that the surface convective cell size could be $\sim 10\%$
of the stellar radius within the uncertainties of parameters measured for
Betelgeuse (Dolan et al., 2016; Joyce et al., 2020).

The convective velocity $v_{c}$ in the RSG envelope can be estimated by
equating the energy flux to the convective flux:

$v_{c}\sim\sqrt[3]{\frac{\sigma_{\mathrm{SB}}T_{\mathrm{eff}}^{4}}{\rho_{s}}}\sim
21\,\mathrm{km\,s^{-1}}\\!\left(\frac{T_{\mathrm{eff}}}{3600\,\mathrm{K}}\right)^{\frac{4}{3}}\rho_{s,-9}^{-\frac{1}{3}}\,,$
(B2)

where $\sigma_{\mathrm{SB}}$ is the Stefan-Boltzmann constant and
$\rho_{s,-9}$ is the surface density in units of
$10^{-9}\,\mathrm{g\,cm^{-3}}$ (representative density taken from the MESA
models; see also Figure 2 of Goldberg et al., 2022). The constants of order
unity are omitted here.

The timescale for the surface convective structure to readjust can be
estimated as

$t_{\mathrm{surf}}\sim\frac{d}{v_{c}}\sim
25\,\mathrm{days}\left(\frac{\alpha}{10}\right)\\!\left(\frac{T_{\mathrm{eff}}}{3600\,\mathrm{K}}\right)^{-\frac{1}{3}}\\!\left(\frac{g}{\mathrm{g\,cm^{-3}}}\right)^{-1}\\!\rho_{s,-9}^{\frac{1}{3}}\,,$
(B3)

typically at the order of weeks to months.

In comparison, the velocity field is predicted to be influenced by the deep
convection, which operates over a larger length scale comparable to the
stellar radius as shown in 3D simulations (Kravchenko et al., 2019). The
timescale for the deep convection can thus be estimated as the convective
turnover timescale

$t_{\mathrm{turnover}}\sim\frac{R}{v_{c}}\sim
0.8\,\mathrm{yrs}\left(\frac{R}{800R_{\odot}}\right)\left(\frac{T_{\mathrm{eff}}}{3600\,\mathrm{K}}\right)^{-\frac{4}{3}}\rho_{s,-9}^{\frac{1}{3}}\,.$
(B4)

##  Appendix C Post-processing package: methods, tests, and limitations

###  C.1 Equilibrium chemistry with FASTCHEM2

Figure 5:  Comparison of number densities from interpolated results of pre-
calculated FASTCHEM2 table (squares) against those from MARCS (curves). The
colors indicate the number densities of atomic hydrogen, electrons, CO
molecules, and SiO molecules, respectively.  Figure 6:  2D middle slice of SiO
number density in the simulated RSG envelope. The green circle is the radio
photosphere of the star. From the same snapshot of the simulation shown in
Figure 2.

The number densities of different species in our simulations are obtained from
FASTCHEM2 assuming chemical equilibrium. On-the-fly calculations are time-
consuming given the large number of grid points and snapshots in the
simulations. Therefore, we pre-calculate the number density tables using
FASTCHEM2 that cover the parameter space of
$2.5\leq\log_{10}(T/[\mathrm{K}])\leq 5.5,\ -30\leq\log_{10}({\cal R})\leq-9$
where ${\cal R}=(P/[\mathrm{dyn\,cm^{-2}}])/(T/[\mathrm{K}])^{4}$. Here, $T$
is the gas temperature and $P$ is the gas pressure. The number densities in
our simulations are interpolated values from these tables. To verify that the
interpolated results from the pre-calculated FASTCHEM2 table are suitable for
RSG simulations, we compare the interpolated values to the values given by a
MARCS RSG model in Figure 5. The MARCS code was developed to create 1D
atmospheric models for evolved stars (Gustafsson et al., 2008), which has been
widely used to produce spectra to obtain basic parameters of RSGs from
spectroscopy (e.g. Levesque et al., 2005). The interpolated FASTCHEM2 values
closely agree with the MARCS values.

In Figure 6, we show an example of the SiO number density distribution in our
simulations by taking a 2D slice along the middle cross-section. The light
green circle indicates the radio photosphere determined by the continuum
intensity map (see Section 2.2). In our simulations, molecules are dissociated
in the inner envelope due to high temperature. The most abundant SiO molecules
are found near the optical photosphere due to recombination. The presence of
shocks outside of the photosphere can be seen from the SiO number density
slice.

###  C.2 Continuum and line radiative transfer

In our post-processing, we numerically integrate the time-independent
radiative transfer equation along the ray neglecting scattering. A distance
$x$ can be defined from the origin of the ray. The optical depth is

$\tau_{\nu}(\hat{\textbf{n}},\textbf{r})\equiv\int_{0}^{x(\hat{\textbf{n}},\textbf{r})}\chi_{\nu}\mathrm{d}x^{\prime}\,.$
(C1)

Here, $\hat{\textbf{n}}$ is the unit vector along the ray direction, and r is
the vector of position. The intensity can then be calculated as

$I_{\nu}(x)=I_{\nu}(0)e^{-\tau_{\nu}}+\int_{e^{-\tau_{\nu}}}^{1}S_{\nu}(\tau_{\nu}-\tau_{\nu}^{\prime})\mathrm{d}(e^{-\tau_{\nu}^{\prime}})\,,$
(C2)

where $S_{\nu}$ is the source function defined by
$S_{\nu}\equiv\eta_{\nu}/\chi_{\nu}$. $I_{\nu},\ \eta_{\nu},\ \chi_{\nu}$ are
the intensity, emissivity, and opacity at frequency $\nu$. The intensity
$I_{\nu}(0)$ at the origin is taken to be zero.

For the wavelength range studied here, the emissivity and opacity are
contributed by continuum emission and molecular lines, namely
$\eta_{\nu}=\eta_{\nu}^{\mathrm{con}}+\eta_{\nu}^{\mathrm{line}},\
\chi_{\nu}=\chi_{\nu}^{\mathrm{con}}+\chi_{\nu}^{\mathrm{line}}$. We assume
that $\eta_{\nu}^{\mathrm{con}}=\chi_{\nu}^{\mathrm{con}}B_{\nu}$, where
$B_{\nu}$ is the Planck function. The continuum opacity
$\chi_{\nu}^{\mathrm{con}}$, dominated by $\mathrm{H^{-}}$ free-free
transitions, is taken from the analytical fits in Harper et al. (2001)
specifically made for radio observations of Betelgeuse. Relevant quantities
for $\mathrm{{}^{28}SiO}$ ($\nu=2,J=8-7$) lines are obtained from the ExoMol
database222https://www.exomol.com (Tennyson & Yurchenko, 2012; Tennyson et
al., 2016) using the results of Yurchenko et al. (2022), and further computed
assuming local thermal equilibrium (LTE). The line emissivity and opacity
contributed by the transition from level i to level j are computed following
standard assumptions (e.g. De Ceuster et al., 2020a):

$\displaystyle\eta_{\nu}^{\mathrm{line},ij}=\frac{h\nu}{4\pi}n_{i}A_{ij}\phi_{\nu}^{ij}\,,$
(C3)
$\displaystyle\chi_{\nu}^{\mathrm{line},ij}=\frac{h\nu}{4\pi}(n_{j}B_{ji}-n_{i}B_{ij})\phi_{\nu}^{ij}\,,$
(C4)

where $h$ is the Planck constant, $n_{i}$ is the number density of the species
at level i, and $A_{ij},\ B_{ij}$ are Einstein coefficients. The line profile
function is assumed to be a Gaussian with a line width contributed by thermal
broadening

$\displaystyle\phi_{\nu}^{ij}=\frac{1}{\delta\nu_{ij}\sqrt{\pi}}\exp{\left[-\left(\frac{\nu-\tilde{\nu}_{ij}}{\delta\nu_{ij}}\right)^{2}\right]}\,,$
(C5)
$\displaystyle\delta\nu_{ij}=\frac{|v_{\mathrm{therm}}|}{c}\tilde{\nu}_{ij}\,,$
(C6)

where the Doppler-shifted frequency $\tilde{\nu}_{ij}$ and the thermal
velocity $v_{\mathrm{therm}}$ are

$\displaystyle\tilde{\nu}_{ij}=\nu_{ij}(1-\frac{\hat{\textbf{n}}\cdot\textbf{v}}{c})\,,$
(C7) $\displaystyle v_{\mathrm{therm}}=\sqrt{\frac{2k_{B}T}{m}}\,.$ (C8)

Here, $\nu_{ij}$ is the rest frequency for the transition, v is the local
velocity, $c$ is the speed of light, $T$ is the gas temperature, and $m$ is
the particle mass of the species. The Doppler shift follows the convention
that velocity towards the observer is blue-shifted and has negative values.

We have tested our algorithm against the line transfer code MAGRITTE (De
Ceuster et al., 2020a, b, 2022) by applying both of them to the CO5BOLD
simulation snapshot. The results agree with each other. However, MAGRITTE
needs to construct the Voronoi grid. In comparison, our method runs faster
with CO5BOLD snapshots which are computed on a Cartesian grid, and does not
suffer from extra interpolation errors. Additionally, MAGRITTE has not yet
incorporated continuum opacity or interface to the ExoMol database, both of
which are important for this study.

###  C.3 Gaussian line fit and rotation measurement

Figure 7:  Examples of emission line profiles from synthetic ALMA images. Four
small subplots show the simulated line profiles (blue) and fitted Gaussian
(orange) in four pixels (shaded in black) of the synthetic ALMA image. In each
pixel, the mean value (orange dashed vertical line) of the fitted Gaussian is
taken to be the radial velocity from the Doppler shift. From the same snapshot
of the simulation shown in Figure 2.

An example of our Gaussian line fit procedure is illustrated in Figure 7.
Green arrows connect the line profiles to their corresponding pixels
(indicated in black). To obtain the radial velocity map, we fit a Gaussian
(orange solid line) to the line profile (blue solid line) in each pixel of the
image, and find the mean value of the Gaussian (orange dotted line).

Following Section 3.5 in K18, we subtract the inferred systematic velocity
from the radial velocity map, and fit a projected radial velocity map of a
rigidly rotating sphere to the synthesized map. Namely, we assume that the
$v_{\mathrm{sys}}$-subtracted radial velocity map has the velocity
distribution of $v(x,y)=v\sin
i\times(x\cos\theta+y\sin\theta)/R_{\mathrm{star}}$, where the fitted
parameters are the position angle $\theta$ and the inferred $v\sin i$, and
$x,y$ are the coordinates of the 2D radial velocity map.

Instead of using the radio photospheric radius $R_{\mathrm{star}}$ for the
fit, we suggest it is more consistent with the physical scenario to use the
radius $R_{\mathrm{shell}}\approx 1.3R_{\mathrm{star}}$ of the molecular shell
probed by ALMA. This would increase the measured $v\sin i$ _at the molecular
shell_ by a factor of $1.3$, i.e., $v\sin i_{\mathrm{shell}}\approx
7\,\mathrm{km\,s^{-1}}$ at the molecular shell. Assuming that the molecular
shell co-rotates with the radio photosphere, this translates to a rotational
velocity at the radio photosphere a factor of $1.3$ less than the value
measured at the molecular shell, i.e., $v\sin i_{\mathrm{radio}}\approx
5.5\,\mathrm{km\,s^{-1}}$ at the radio photosphere. However, if the molecular
shell and the radio photosphere are not fully coupled, e.g. assuming similar
specific angular momentum between the molecular shell and the radio
photosphere, then the rotational velocity at the radio photosphere would be
another factor of $1.3$ higher, namely $v\sin i_{\mathrm{radio}}\approx
9\,\mathrm{km\,s^{-1}}$ at the radio photosphere. Since the optical
photosphere is yet another factor of $1.3$ smaller than the radio photosphere,
the measurement by K18 would give $v\sin i_{\mathrm{optical}}\approx
4\,\mathrm{km\,s^{-1}}$ at the optical photosphere assuming constant angular
frequency (co-rotation) outside the star, or $v\sin
i_{\mathrm{optical}}\approx 12\,\mathrm{km\,s^{-1}}$ assuming constant
specific angular momentum outside the star. In this work, we closely follow
the procedure done by K18 for direct comparison, and thus use
$R_{\mathrm{star}}$ for the rotation fit, with the underlying assumption that
the molecular shell co-rotates with the radio photosphere.

Table 2: Parameters of the CO5BOLD simulations taken from Ahmad et al. (2023), compared to the observed values for Betelgeuse. The stellar parameters for both simulations and observations are listed as surface gravity $\log_{10}g$, effective temperature $T_{\mathrm{eff}}$, luminosity $\log_{10}L_{\star}$, and radius $R_{\star}$. All these surface quantities in simulations are taken at the local minimum point of specific entropy averaged over spherical shells and time, as described in Ahmad et al. (2023). We also list the mass $M_{\mathrm{pot}}$ used for the gravitational potential and envelope mass $M_{\mathrm{env}}$ for the simulations. The convective turnover timescale $t_{\mathrm{turnover}}$ is estimated analytically as $R_{\star}/v_{c}$, where the convective velocity $v_{c}$ are calculated using the equations in Appendix B. The number of grid cells, simulation box size $x_{\mathrm{box}}$, number of snapshots, and the physical time $t_{\mathrm{snap}}$ covered by the extracted snapshots are listed in sequence.  |  | $\log_{10}g$ | $T_{\mathrm{eff}}$ | $\log_{10}L_{\star}$ | $R_{\star}$ | $M_{\mathrm{pot}}$ | $M_{\mathrm{env}}$ | $t_{\mathrm{turnover}}$ | Grid | $x_{\mathrm{box}}$ | Snapshot | $t_{\mathrm{snap}}$
---|---|---|---|---|---|---|---|---|---|---|---|---
|  |  [$\mathrm{cm\,s^{-2}}$] | [K] | $[L_{\odot}]$ | $[R_{\odot}]$ | $[M_{\odot}]$ | $[M_{\odot}]$ | [yr] | [#] | $[R_{\odot}]$ | [#] | [yr]
Sim. | default | $-$0.419 | 3366 | 4.62 | 597 | 5.0 | 0.5 | $\sim 0.7$ | 3153 | 1626 | 80 | 5.1
(st35gm04n37)
B | $-$0.246 | 3620 | 4.95 | 759 | 12.0 | 3.0 | $\sim 0.8$ | 6373 | 2093 | 100 | 6.3
(st35gm03n020)
Obs. | Betelgeuse |  $-0.5$aaLobel & Dupree (2000) |  $3500^{+200}_{-200}$bbDolan et al. (2016) (chosen value motivated by multiple sources) |  $5.10^{+0.22}_{-0.22}$ccHarper et al. (2008) |  $764^{+116}_{-62}$ddJoyce et al. (2020) | - | - | - | - | - | - | -
Figure 8:  Continuum intensity map and radial velocity map of best-fit
simulation B snapshot, compared with ALMA observations of Betelgeuse. Physical
quantities marked in white are defined the same way as in Table 1 using the
method in Section 2.3. Credit for the right two panels: Kervella et al., A&A,
609, A67, 2018, reproduced with permission © ESO.

###  C.4 Limitations

We scale the simulation such that it has an inferred radio photospheric radius
$R_{\mathrm{star}}$ of $\sim 30$ mas. However, this corresponds to a distance
of 130 pc, consistent with the original _Hipparcos_ parallax
($131^{+35}_{-23}$; ESA, 1997; Perryman et al., 1997) and its revised value
($152^{+22}_{-17}$; Van Leeuwen, 2007) but much closer than the distance
derived from _Hipparcos_ data combined with radio observations
($197^{+55}_{-35}$; Harper et al., 2008) and its revised value
($222^{+48}_{-34}$; Harper et al., 2017). This is partly because the default
simulation used in this work has a smaller radius than Betelgeuse, and
therefore needs to be placed closer to get a similar parallax. Another reason
is rooted in the limitations of the 3D models, which have been shown to be
less extended than actual RSGs (Arroyo-Torres et al., 2015; Climent et al.,
2020; Chiavassa et al., 2022). This could be due to the fact that the
radiation pressure is not included in the simulations (Freytag et al., 2012).

The simulated continuum intensity is about half of the observed value (see
Figure 2), and the line amplitude is an order of magnitude lower than observed
(see Figure 7 compared to Figure 1 and 9 in K18). This could result from the
small atmospheric extension of the simulations, or not including non-
equilibrium chemistry, non-LTE populations, dust emission, maser, or
scattering in the post-processing. However, as we are interested in the radial
velocity shifts instead of the molecular abundances in this study, the effects
of these missing physics are expected to be limited. This is supported by our
simulations, which show no significant variations in the velocity field from
the surface to the molecular layer.

A more uncertain aspect is whether the convective structure in our simulations
can represent actual RSG convection. The convective length scales are subject
to uncertainties of stellar parameters of Betelgeuse itself, as well as the
resolution of the simulations (Chiavassa et al., 2011a). Aside from that, core
convection simulations suggest that large-scale convection is dominated by
dipolar flows (see, e.g. Lecoanet & Edelmann, 2023, for a review), which are
specifically damped in our simulations by the drag forces in the central
region (Freytag et al., 2012). However, we expect that an enhancement in the
dipolar flows would only make our argument stronger.

Figure 9:  Probability distribution of apparent $v\sin i$, similar to Figure 3
but inferred from simulation B.

##  Appendix D Additional simulations

###  D.1 Impact of the chosen stellar parameters

The parameters of two sets of 3D simulations used in this study are shown in
Table 2. The default simulation is discussed in the main text (Figure 2 and
Figure 3), while an alternative simulation with smaller surface granules is
presented in Figure 8 and Figure 9. Due to the small granule size compared to
the ALMA beam size, the radial velocity map displays small-scale features
(Figure 8) and therefore is less likely to be mistaken for rotation (Figure
9). The change in the intensity map, however, is less obvious. In particular,
its radial velocity map shows a velocity magnitude similar to the default
simulation, but with more cell-like structures. This indicates that the
magnitude of velocity in the radial velocity map may be only determined by the
maximum convective velocity and the relative beam size, and not sensitive to
the cell size. The morphology of the radial velocity map, however, is
sensitive to the cell size, and consequently to the surface gravity.

###  D.2 Impact of the beam size

To further explore the importance of the beam size relative to the surface
granule size, we vary the beam size in the mock observations for the default
simulation. In Figure 10, we show that, as the spatial resolution of the
interferometer increases, more small-scale features can be revealed in both
the intensity map and the radial velocity map. As shown in the third column (c
& g), decreasing the beam size by a factor of two would already help to
marginally resolve the convective structure, both in the radial velocity map
and the intensity map. On the contrary, increasing the beam size by a factor
of two (a & e) smears out any asymmetrical feature in the intensity map, and
the radial velocity map is completely dominated by a dipolar pattern.
Strikingly, as the spatial resolution of the interferometer increases (left to
right), the measured radio photospheric radius monotonically decreases from
$\sim 30$ mas to $\sim 25$ mas, whilst the peak magnitude of radial velocity
monotonically increases from $\sim 6$ to $\sim 30\,\mathrm{km\,s^{-1}}$. As a
result of increasing radial velocity magnitude, the fitted $v\sin i$,
$\chi^{2}_{\mathrm{red}}$, and $\sigma_{\mathrm{res}}$ all increase with
higher resolution, among which $\chi^{2}_{\mathrm{red}}$ increases by three
orders of magnitude. If future observations can be done with higher
resolution, these are the clear signatures to look for to support our
hypothesis.

Figure 10:  Effects of spatial resolution on the observed intensity map and
radial velocity map in mock observations. From the left to the right, spatial
resolution increases with decreasing beam size (indicated by grey circles).
The second column (b & f) uses the standard resolution of FWHM $=18$ mas
adopted in ALMA observations of Betelgeuse (K18), the same as shown in Figure
2. In comparison, the beam size of the first column (a & e) is increased by a
factor of two (36 mas), whereas the beam size of the third column (c & g) is
reduced by a factor of two (9 mas). The last column (d & h) is the same
original simulation image shown in Figure 2. Physical quantities marked in
white are defined the same way as in Table 1 using the method in Section 2.3.

## References

  * Agúndez et al. (2020) Agúndez, M., Martínez, J. I., de Andres, P. L., Cernicharo, J., & Martín-Gago, J. A. 2020, A&A, 637, A59, doi: 10.1051/0004-6361/202037496
  * Ahmad et al. (2023) Ahmad, A., Freytag, B., & Höfner, S. 2023, A&A, 669, A49, doi: 10.1051/0004-6361/202244555
  * Andriantsaralaza et al. (2021) Andriantsaralaza, M., Ramstedt, S., Vlemmings, W. H. T., et al. 2021, A&A, 653, A53, doi: 10.1051/0004-6361/202140952
  * Antoni & Quataert (2022) Antoni, A., & Quataert, E. 2022, MNRAS, 511, 176, doi: 10.1093/mnras/stab3776
  * Ariste et al. (2022) Ariste, A. L., Georgiev, S., Mathias, P., et al. 2022, A&A, 661, A91, doi: 10.1051/0004-6361/202142271
  * Ariste et al. (2018) Ariste, A. L., Mathias, P., Tessore, B., et al. 2018, A&A, 620, A199, doi: 10.1051/0004-6361/201834178
  * Arroyo-Torres et al. (2015) Arroyo-Torres, B., Wittkowski, M., Chiavassa, A., et al. 2015, A&A, 575, A50, doi: 10.1051/0004-6361/201425212
  * Asaki et al. (2023) Asaki, Y., Maud, L. T., Francke, H., et al. 2023, ALMA High-frequency Long Baseline Campaign in 2021: Highest Angular Resolution Submillimeter Wave Images for the Carbon-rich Star R Lep, doi: 10.48550/arXiv.2310.09664
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Barnbaum et al. (1995) Barnbaum, C., Morris, M., & Kahane, C. 1995, ApJ, 450, 862, doi: 10.1086/176190
  * Brott et al. (2011) Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115, doi: 10.1051/0004-6361/201016113
  * Brun & Palacios (2009) Brun, A. S., & Palacios, A. 2009, ApJ, 702, 1078, doi: 10.1088/0004-637X/702/2/1078
  * Brun & Toomre (2002) Brun, A. S., & Toomre, J. 2002, ApJ, 570, 865, doi: 10.1086/339228
  * Brunner et al. (2019) Brunner, M., Mecina, M., Maercker, M., et al. 2019, A&A, 621, A50, doi: 10.1051/0004-6361/201833652
  * Carlberg et al. (2012) Carlberg, J. K., Cunha, K., Smith, V. V., & Majewski, S. R. 2012, ApJ, 757, 109, doi: 10.1088/0004-637X/757/2/109
  * Carpenter et al. (2020) Carpenter, J., Iono, D., Kemper, F., & Wootten, A. 2020, The ALMA Development Program: Roadmap to 2030, doi: 10.48550/arXiv.2001.11076
  * Chatzopoulos et al. (2020) Chatzopoulos, E., Frank, J., Marcello, D. C., & Clayton, G. C. 2020, ApJ, 896, 50, doi: 10.3847/1538-4357/ab91bb
  * Chiavassa et al. (2011a) Chiavassa, A., Freytag, B., Masseron, T., & Plez, B. 2011a, A&A, 535, A22, doi: 10.1051/0004-6361/201117463
  * Chiavassa et al. (2010) Chiavassa, A., Haubois, X., Young, J. S., et al. 2010, A&A, 515, A12, doi: 10.1051/0004-6361/200913907
  * Chiavassa et al. (2024) Chiavassa, A., Kravchenko, K., & Goldberg, J. A. 2024, Living Reviews in Computational Astrophysics, submitted
  * Chiavassa et al. (2009) Chiavassa, A., Plez, B., Josselin, E., & Freytag, B. 2009, A&A, 506, 1351, doi: 10.1051/0004-6361/200911780
  * Chiavassa et al. (2011b) Chiavassa, A., Pasquato, E., Jorissen, A., et al. 2011b, A&A, 528, A120, doi: 10.1051/0004-6361/201015768
  * Chiavassa et al. (2022) Chiavassa, A., Kravchenko, K., Montargès, M., et al. 2022, A&A, 658, A185, doi: 10.1051/0004-6361/202142514
  * Climent et al. (2020) Climent, J. B., Wittkowski, M., Chiavassa, A., et al. 2020, A&A, 635, A160, doi: 10.1051/0004-6361/201936734
  * De Ceuster et al. (2020a) De Ceuster, F., Homan, W., Yates, J., et al. 2020a, MNRAS, 492, 1812, doi: 10.1093/mnras/stz3557
  * De Ceuster et al. (2020b) De Ceuster, F., Bolte, J., Homan, W., et al. 2020b, MNRAS, 499, 5194, doi: 10.1093/mnras/staa3199
  * De Ceuster et al. (2022) De Ceuster, F., Ceulemans, T., Srivastava, A., et al. 2022, JOSS, 7, 3905, doi: 10.21105/joss.03905
  * de Jager et al. (1988) de Jager, C., Nieuwenhuijzen, H., & van der Hucht, K. A. 1988, A&AS, 72, 259.  https://ui.adsabs.harvard.edu/abs/1988A%26AS...72..259D/abstract
  * de Mink et al. (2014) de Mink, S. E., Sana, H., Langer, N., Izzard, R. G., & Schneider, F. R. N. 2014, ApJ, 782, 7, doi: 10.1088/0004-637X/782/1/7
  * Decin et al. (2023) Decin, L., Richards, A. M. S., Marchant, P., & Sana, H. 2023, ALMA detection of CO rotational line emission in red supergiant stars of the massive young star cluster RSGC1 – Determination of a new mass-loss rate prescription for red supergiants, doi: 10.48550/arXiv.2303.09385
  * Decin et al. (2012) Decin, L., Cox, N. L. J., Royer, P., et al. 2012, A&A, 548, A113, doi: 10.1051/0004-6361/201219792
  * Decin et al. (2020) Decin, L., Montargès, M., Richards, A. M. S., et al. 2020, Science, 369, 1497, doi: 10.1126/science.abb1229
  * Doan et al. (2020) Doan, L., Ramstedt, S., Vlemmings, W. H. T., et al. 2020, A&A, 633, A13, doi: 10.1051/0004-6361/201935245
  * Dolan et al. (2016) Dolan, M. M., Mathews, G. J., Lam, D. D., et al. 2016, ApJ, 819, 7, doi: 10.3847/0004-637X/819/1/7
  * Drevon et al. (2024) Drevon, J., Millour, F., Cruzalèbes, P., et al. 2024, MNRAS, 527, L88, doi: 10.1093/mnrasl/slad138
  * ESA (1997) ESA. 1997, ESA SP Series, Vol. 1200, The HIPPARCOS and TYCHO catalogues. Astrometric and photometric star catalogues derived from the ESA HIPPARCOS Space Astrometry Mission (Noordwijk, Netherlands: ESA Publications Division).  https://ui.adsabs.harvard.edu/abs/1997ESASP1200.....E
  * Freytag et al. (2012) Freytag, B., Steffen, M., Ludwig, H. G., et al. 2012, JCoPh, 231, 919, doi: 10.1016/j.jcp.2011.09.026
  * Gaulme et al. (2020) Gaulme, P., Jackiewicz, J., Spada, F., et al. 2020, A&A, 639, A63, doi: 10.1051/0004-6361/202037781
  * Gilliland & Dupree (1996) Gilliland, R. L., & Dupree, A. K. 1996, ApJ, 463, L29, doi: 10.1086/310043
  * Goldberg et al. (2022) Goldberg, J. A., Jiang, Y.-F., & Bildsten, L. 2022, ApJ, 929, 156, doi: 10.3847/1538-4357/ac5ab3
  * Gottlieb et al. (2022) Gottlieb, C. A., Decin, L., Richards, A. M. S., et al. 2022, A&A, 660, A94, doi: 10.1051/0004-6361/202140431
  * Guinan et al. (2020) Guinan, E., Wasatonic, R., Calderwood, T., & Carona, D. 2020, The Astronomer’s Telegram, 13512, 1.  https://ui.adsabs.harvard.edu/abs/2020ATel13512....1G
  * Guinan et al. (2019) Guinan, E. F., Wasatonic, R. J., & Calderwood, T. J. 2019, The Astronomer’s Telegram, 13341, 1.  https://ui.adsabs.harvard.edu/abs/2019ATel13341....1G
  * Gustafsson et al. (2008) Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951, doi: 10.1051/0004-6361:200809724
  * Harper & Brown (2006) Harper, G. M., & Brown, A. 2006, ApJ, 646, 1179, doi: 10.1086/505073
  * Harper et al. (2008) Harper, G. M., Brown, A., & Guinan, E. F. 2008, AJ, 135, 1430, doi: 10.1088/0004-6256/135/4/1430
  * Harper et al. (2017) Harper, G. M., Brown, A., Guinan, E. F., et al. 2017, ApJ, 154, 11, doi: 10.3847/1538-3881/aa6ff9
  * Harper et al. (2001) Harper, G. M., Brown, A., & Lim, J. 2001, ApJ, 551, 1073, doi: 10.1086/320215
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Haubois et al. (2009) Haubois, X., Perrin, G., Lacour, S., et al. 2009, A&A, 508, 923, doi: 10.1051/0004-6361/200912927
  * Heger et al. (2000) Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368, doi: 10.1086/308158
  * Heger et al. (2005) Heger, A., Woosley, S. E., & Spruit, H. C. 2005, ApJ, 626, 350, doi: 10.1086/429868
  * Herwig (2000) Herwig, F. 2000, A&A, 360, 952, doi: 10.48550/arXiv.astro-ph/0007139
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jadlovský et al. (2023) Jadlovský, D., Krtička, J., Paunzen, E., & Štefl, V. 2023, New Astronomy, 99, 101962, doi: 10.1016/j.newast.2022.101962
  * Jermyn et al. (2023) Jermyn, A. S., Bauer, E. B., Schwab, J., et al. 2023, ApJS, 265, 15, doi: 10.3847/1538-4365/acae8d
  * Josselin & Plez (2007) Josselin, E., & Plez, B. 2007, A&A, 469, 671, doi: 10.1051/0004-6361:20066353
  * Joyce et al. (2020) Joyce, M., Leung, S.-C., Molnár, L., et al. 2020, ApJ, 902, 63, doi: 10.3847/1538-4357/abb8db
  * Kee et al. (2021) Kee, N. D., Sundqvist, J. O., Decin, L., Koter, A. d., & Sana, H. 2021, A&A, 646, A180, doi: 10.1051/0004-6361/202039224
  * Kervella et al. (2018) Kervella, P., Decin, L., Richards, A. M. S., et al. 2018, A&A, 609, A67, doi: 10.1051/0004-6361/201731761
  * Kitzmann et al. (2023) Kitzmann, D., Hoeijmakers, H. J., Grimm, S. L., et al. 2023, A&A, 669, A113, doi: 10.1051/0004-6361/202142969
  * Kitzmann et al. (2018) Kitzmann, D., Heng, K., Rimmer, P. B., et al. 2018, ApJ, 863, 183, doi: 10.3847/1538-4357/aace5a
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas (IOS Press), doi: 10.3233/978-1-61499-649-1-87
  * Kochanek et al. (2014) Kochanek, C. S., Adams, S. M., & Belczynski, K. 2014, MNRAS, 443, 1319, doi: 10.1093/mnras/stu1226
  * Kravchenko et al. (2019) Kravchenko, K., Chiavassa, A., Eck, S. V., et al. 2019, A&A, 632, A28, doi: 10.1051/0004-6361/201935809
  * Kravchenko et al. (2021) Kravchenko, K., Jorissen, A., Van Eck, S., et al. 2021, A&A, 650, L17, doi: 10.1051/0004-6361/202039801
  * Langer et al. (1985) Langer, N., El Eid, M. F., & Fricke, K. J. 1985, A&A, 145, 179.  https://ui.adsabs.harvard.edu/abs/1985A&A...145..179L
  * Langer et al. (1983) Langer, N., Fricke, K. J., & Sugimoto, D. 1983, A&A, 126, 207.  https://ui.adsabs.harvard.edu/abs/1983A&A...126..207L/abstract
  * Lau et al. (2022) Lau, M. Y. M., Cantiello, M., Jermyn, A. S., et al. 2022, Hot Jupiter engulfment by a red giant in 3D hydrodynamics, doi: 10.48550/arXiv.2210.15848
  * Lecoanet & Edelmann (2023) Lecoanet, D., & Edelmann, P. V. F. 2023, Galaxies, 11, 89, doi: 10.3390/galaxies11040089
  * Ledoux (1947) Ledoux, P. 1947, ApJ, 105, 305, doi: 10.1086/144905
  * Levesque et al. (2005) Levesque, E. M., Massey, P., Olsen, K. A. G., et al. 2005, ApJ, 628, 973, doi: 10.1086/430901
  * Lobel & Dupree (2000) Lobel, A., & Dupree, A. K. 2000, ApJ, 545, 454, doi: 10.1086/317784
  * Lobel & Dupree (2001) —. 2001, ApJ, 558, 815, doi: 10.1086/322284
  * Maeder & Meynet (2000) Maeder, A., & Meynet, G. 2000, ARA&A, 38, 143, doi: 10.1146/annurev.astro.38.1.143
  * Mestel (1968) Mestel, L. 1968, MNRAS, 138, 359, doi: 10.1093/mnras/138.3.359
  * Michelson & Pease (1921) Michelson, A. A., & Pease, F. G. 1921, ApJ, 53, 249, doi: 10.1086/142603
  * Montargès et al. (2014) Montargès, M., Kervella, P., Perrin, G., et al. 2014, A&A, 572, A17, doi: 10.1051/0004-6361/201423538
  * Montargès et al. (2018) Montargès, M., Norris, R., Chiavassa, A., et al. 2018, A&A, 614, A12, doi: 10.1051/0004-6361/201731471
  * Montargès et al. (2016) Montargès, M., Kervella, P., Perrin, G., et al. 2016, A&A, 588, A130, doi: 10.1051/0004-6361/201527028
  * Montargès et al. (2023) Montargès, M., Cannon, E., de Koter, A., et al. 2023, A&A, 671, A96, doi: 10.1051/0004-6361/202245398
  * Nhung et al. (2023) Nhung, P. T., Hoai, D. T., Darriulat, P., et al. 2023, RAA, 23, 015004, doi: 10.1088/1674-4527/aca06f
  * Nhung et al. (2021) Nhung, P. T., Hoai, D. T., Tuan-Anh, P., et al. 2021, MNRAS, 504, 2687, doi: 10.1093/mnras/stab954
  * Norris et al. (2021) Norris, R. P., Baron, F. R., Monnier, J. D., et al. 2021, ApJ, 919, 124, doi: 10.3847/1538-4357/ac0c7e
  * O’Gorman et al. (2017) O’Gorman, E., Kervella, P., Harper, G. M., et al. 2017, A&A, 602, L10, doi: 10.1051/0004-6361/201731171
  * Ohnaka et al. (2013) Ohnaka, K., Hofmann, K. H., Schertl, D., et al. 2013, A&A, 555, A24, doi: 10.1051/0004-6361/201321063
  * Ohnaka et al. (2017) Ohnaka, K., Weigelt, G., & Hofmann, K. H. 2017, Nature, 548, 310, doi: 10.1038/nature23445
  * Ohnaka et al. (2019) Ohnaka, K., Weigelt, G., & Hofmann, K.-H. 2019, ApJ, 883, 89, doi: 10.3847/1538-4357/ab3d2a
  * Ohnaka et al. (2009) Ohnaka, K., Hofmann, K.-H., Benisty, M., et al. 2009, A&A, 503, 183, doi: 10.1051/0004-6361/200912247
  * Ohnaka et al. (2011) Ohnaka, K., Weigelt, G., Millour, F., et al. 2011, A&A, 529, A163, doi: 10.1051/0004-6361/201016279
  * Paladini et al. (2018) Paladini, C., Baron, F., Jorissen, A., et al. 2018, Nature, 553, 310, doi: 10.1038/nature25001
  * Patton et al. (2023) Patton, R. A., Pinsonneault, M. H., Cao, L., et al. 2023, Spectroscopic identification of rapidly rotating red giant stars in APOKASC-3 and APOGEE DR16, doi: 10.48550/arXiv.2303.08151
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3, doi: 10.1088/0067-0049/192/1/3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4, doi: 10.1088/0067-0049/208/1/4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15, doi: 10.1088/0067-0049/220/1/15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34, doi: 10.3847/1538-4365/aaa5a8
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10, doi: 10.3847/1538-4365/ab2241
  * Perryman et al. (1997) Perryman, M. A. C., Lindegren, L., Kovalevsky, J., et al. 1997, A&A, 323, L49.  https://ui.adsabs.harvard.edu/abs/1997A&A...323L..49P
  * Privitera et al. (2016a) Privitera, G., Meynet, G., Eggenberger, P., et al. 2016a, A&A, 591, A45, doi: 10.1051/0004-6361/201528044
  * Privitera et al. (2016b) —. 2016b, A&A, 593, A128, doi: 10.1051/0004-6361/201628758
  * Ramstedt et al. (2017) Ramstedt, S., Mohamed, S., Vlemmings, W. H. T., et al. 2017, A&A, 605, A126, doi: 10.1051/0004-6361/201730934
  * Ramstedt et al. (2020) Ramstedt, S., Vlemmings, W. H. T., Doan, L., et al. 2020, A&A, 640, A133, doi: 10.1051/0004-6361/201936874
  * Rui & Fuller (2021) Rui, N. Z., & Fuller, J. 2021, MNRAS, 508, 1618, doi: 10.1093/mnras/stab2528
  * Sana et al. (2012) Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444, doi: 10.1126/science.1223344
  * Schwarzschild (1975) Schwarzschild, M. 1975, ApJ, 195, 137, doi: 10.1086/153313
  * Shiber et al. (2023) Shiber, S., Chatzopoulos, E., Munson, B., & Frank, J. 2023, Betelgeuse as a Merger of a Massive Star with a Companion.  https://ui.adsabs.harvard.edu/abs/2023arXiv231014603S
  * Smith (2014) Smith, N. 2014, ARA&A, 52, 487, doi: 10.1146/annurev-astro-081913-040025
  * Spruit (2002) Spruit, H. C. 2002, A&A, 381, 923, doi: 10.1051/0004-6361:20011465
  * Stock et al. (2022) Stock, J. W., Kitzmann, D., & Patzer, A. B. C. 2022, MNRAS, 517, 4070, doi: 10.1093/mnras/stac2623
  * Stock et al. (2018) Stock, J. W., Kitzmann, D., Patzer, A. B. C., & Sedlmayr, E. 2018, MNRAS, 479, 865, doi: 10.1093/mnras/sty1531
  * Sullivan et al. (2020) Sullivan, J. M., Nance, S., & Wheeler, J. C. 2020, ApJ, 905, 128, doi: 10.3847/1538-4357/abc3c9
  * Tennyson & Yurchenko (2012) Tennyson, J., & Yurchenko, S. N. 2012, MNRAS, 425, 21, doi: 10.1111/j.1365-2966.2012.21440.x
  * Tennyson et al. (2016) Tennyson, J., Yurchenko, S. N., Al-Refaie, A. F., et al. 2016, JMoSp, 327, 73, doi: 10.1016/j.jms.2016.05.002
  * Tremblay et al. (2013) Tremblay, P. E., Ludwig, H. G., Freytag, B., Steffen, M., & Caffau, E. 2013, A&A, 557, A7, doi: 10.1051/0004-6361/201321878
  * Uitenbroek et al. (1998) Uitenbroek, H., Dupree, A. K., & Gilliland, R. L. 1998, AJ, 116, 2501, doi: 10.1086/300596
  * Van Leeuwen (2007) Van Leeuwen, F., ed. 2007, Astrophysics and Space Science Library, Vol. 350, Hipparcos, the New Reduction of the Raw Data (Dordrecht: Springer Netherlands), doi: 10.1007/978-1-4020-6342-8
  * Vink et al. (2001) Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574, doi: 10.1051/0004-6361:20010127
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Vlemmings et al. (2018) Vlemmings, W. H. T., Khouri, T., De Beck, E., et al. 2018, A&A, 613, L4, doi: 10.1051/0004-6361/201832929
  * Wheeler & Chatzopoulos (2023) Wheeler, J. C., & Chatzopoulos, E. 2023, Astronomy & Geophysics, 64, 3.11, doi: 10.1093/astrogeo/atad020
  * Wheeler et al. (2017) Wheeler, J. C., Nance, S., Diaz, M., et al. 2017, MNRAS, 465, 2654, doi: 10.1093/mnras/stw2893
  * Yoon & Cantiello (2010) Yoon, S.-C., & Cantiello, M. 2010, ApJ, 717, L62, doi: 10.1088/2041-8205/717/1/L62
  * Young et al. (2000) Young, J. S., Baldwin, J. E., Boysen, R. C., et al. 2000, MNRAS, 315, 635, doi: 10.1046/j.1365-8711.2000.03438.x
  * Yurchenko et al. (2022) Yurchenko, S. N., Tennyson, J., Syme, A.-M., et al. 2022, MNRAS, 510, 903, doi: 10.1093/mnras/stab3267
  * Zapartas et al. (2019) Zapartas, E., de Mink, S. E., Justham, S., et al. 2019, A&A, 631, A5, doi: 10.1051/0004-6361/201935854


Convolutional Matching Pursuit and Dictionary Training

Arthur Szlam, Koray Kavukcuoglu, and Yann LeCun

##  1 Introduction

One of the most succesful recent signal processing paradigms has been the
sparse coding/dictionary design model [8, 4]. In this model, we try to
represent a given $d\times n$ data matrix $X$ of $n$ points in
$\mathbb{R}^{d}$ written as columns via a solution to the problem

$\\{W_{*},Z_{*}\\}=\\{W_{*}(K,X,q),Z_{*}(K,X,q)\\}$
$=arg\min\limits_{Z\in\mathbb{R}^{K\times n},W\in\mathbb{R}^{d\times
K}}\sum_{k}||Wz_{k}-x_{k}||^{2},\,\,||z_{k}||_{0}\leq q,$ (1.1)

or its $Z$ coordinate convexification

$\\{\tilde{W}_{*},\tilde{Z}_{*}\\}=\\{\tilde{W}_{*}(K,X,\lambda),\tilde{Z}_{*}(K,X,\lambda)\\}$
$=arg\min\limits_{Z\in\mathbb{R}^{K\times n},W\in\mathbb{R}^{d\times
K}}\sum_{k}||Wz_{k}-x_{k}||^{2}+\lambda||z_{k}||_{1}.$ (1.2)

Here, $\\{W,Z\\}$ are the dictionary and the coefficients, respectively, and
$z_{k}$ is the $k$th column of $Z$. $K$, $q$, and $\lambda$ are user selected
parameters controlling the power of the model.

More recently, many models with additional structure have been proposed. For
example, in [9, 2], the dictionary elements are arranged in groups and the
sparsity is on the group level. In [3, 5, 7], the dictionaries are constructed
to be translation invariant. In the former work, the dictionary is constructed
via a non-negative matrix factorization. In the latter two works, the
construction is a convolutional analogue of 1.2 or an $l^{p}$ variant, with
$0<p<1$. In this short note we work with greedy algorithms for solving the
convolutional analogues of 1.1. Specifically, we demonstrate that sparse
coding by matching pursuit and dictionary learning via K-SVD [1] can be used
in the translation invariant setting.

##  2 Matching Pursuit

Matching pursuit [6] is a greedy algorithm for the solution of the sparse
coding problem

$\min_{z}||Wz-x||^{2},$ $||z||_{0}\leq q,$

where the $d\times k$ matrix $W$ is the dictionary, the $k\times 1$ $z$ is the
code, and $x$ is an $d\times 1$ data vector.

  1. 1.

Set $e=x$, and $z$ the $k$-dimensional zero vector.

  2. 2.

Find $j={\rm arg}\max\limits_{i}||W_{i}^{T}e||_{2}^{2}$.

  3. 3.

Set $a=W_{j}^{T}x$.

  4. 4.

Set $e\leftarrow e-aW_{j}$, and $z_{j}=z_{j}+a$.

  5. 5.

Repeat for $q$ steps

Note that with a bit of bookkeeping, it is only necessary to multiply $W$
against $x$ once, instead of $q$ times. This at a cost of an extra $O(K^{2})$
storage: set $e_{r}$ and $a_{r}$ be $e$ and $a$ from the $r$th step above.
Then:

$W^{T}e_{0}=W^{T}x;$ $W^{T}e_{1}=W^{T}x-a_{0}W^{T}W_{j_{0}},$

and so on. If the Gram matrix for $W$ is stored, this is just a lookup.

###  2.1 Convolutional MP

We consider the special case

$\min_{z}||\sum_{j=1}^{k}w_{j}*z_{j}-x||^{2},$ $||\overline{z}||_{0}\leq q,$

where each $w_{j}$ is a filter, and $\overline{z}$ is all of the responses.

Note that the Gram matrix of the “Toeplitz” dictionary consisting of all the
shifts of the $w_{j}$ is usually too big to be used as a lookup table.
However, because of the symmetries of the convolution, it is also unnecessary;
we only need store a $4*h_{f}\times w_{f}\times k^{2}$ array of inner
products, where $h_{f}$ and $w_{f}$ are the dimensions of the filters.

With this additional storage, to run $q$ basis pursuit steps with $k$ filters
on an $h\times w$ image costs the computation of one application of the filter
bank plus $O(kqhw)$ operations.

##  3 Learning the filters

Given a set of $x$, we can learn the filters and the codes simultaneously.
Several methods are available. A simple one is to alternate between updating
the codes and updating the filters, as in K-SVD [1]:

  1. 1.

Initialize $k$ $h_{f}\times w_{f}$ filters $\\{w_{1},...,w_{k}\\}$.

  2. 2.

Solve for $z$ as above.

  3. 3.

For each filter $w_{j}$,

     * •

find all locations in all the data images where $w_{j}$ is activated

     * •

extract the $h_{f}\times w_{f}$ patch $E_{p}$ from the reconstruction via $z$
at each activated point $p$.

     * •

remove the contribution of $w_{j}$ from each $E_{p}$ (i.e. $E_{p}\leftarrow
E_{p}-c_{(p,j)}w_{j}$, where $c_{(p,j)}$ was the activation determined by
$z$).

     * •

update $w_{j}\leftarrow\text{PCA}(E_{p})$

  4. 4.

Repeat from step 2 until fixed number of iterations.

We note that the forward subproblem (finding $Z$ with $W$ fixed) is not
convex, and so the alternation is not guaranteed to decrease the energy or to
converge to even a local minimum. However, in practice, on image and audio
data, this method generates good filters.

##  4 Some experiments

We train filters on three data sets: the AT&T face database, the motorcycles
from a Caltech database, and the VOC PASCAL database. For all the images in
all our experiments, we perform an additive contrast normalization: each image
$x$ is transformed into $x^{\prime}=x-x*b$, where $b$ is a $5\times 5$
averaging box filter. This is very nearly transforming
$x^{\prime}=\nabla^{2}x$, that is, using the discrete Laplacian of the image
instead of the image. Using the Laplacian would correspond to using the energy

$\sum_{x}||\nabla\left(\sum_{j}w_{j}*z_{j}-x\right)||^{2},$

that is, the energy sees the difference between gradients, not intensities.

###  4.1 Faces

The AT&T face database, available at
http://www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase.html is a set of
400 images of 40 individuals. The faces are centered in each image. We resize
each image to $64\times 64$ and contrast normalize. We train $8$ $16\times 16$
filters. After training the filters we find the feature maps of each image in
the database, obtaining a new set of 400 $8$ channel images. We take the
elementwise absolute value of each of the 8 channel images, and then average
pool over $8\times 8$ blocks. We then train a new 16 element dictionary on the
subsampled images. In figure 1 we display the first level filters, and the
second level filters up to shifts of size 8 and sign changes of the first
level filters..

Figure 1: First and second layer filters from faces

Figure 2: a contrast normalized face, and its reconstruction from 40 filter
responses.

###  4.2 Caltech motorcycles

We also train on the motorbikes-side dataset, available at
http://www.vision.caltech.edu/html-files/archive.html which consists of color
images of various motorcycles. The motorcycles are centered in each image. We
convert each image to gray level, resize to $64\times 64$, and contrast
normalize. We train $8$ $16\times 16$ filters. As before, we then train a new
16 element dictionary on the subsampled absolute value rectified responses of
the first level. In figure 3 we display the first level filters, and the
second level filters up to shifts of size 8 and sign changes of the first
level filters..

Figure 3: First and second layer filters from motorcycles

Figure 4: A contrast normalized motorcycle, and its reconstruction from 40
filter responses.

###  4.3 Images from PASCAL VOC

We also show results trained on “unclassified” natural images from the PASCAL
visual object challenge dataset available at
http://pascallin.ecs.soton.ac.uk/challenges/VOC/. We randomly subsample $5000$
grayscaled images by a factor of 1 to 4, and then pick from each image a
$64\times 64$ patch, and then contrast normalize. We train $8$ $8\times 8$
filters. We then train a new $4\times 4$ 64 element dictionary on the
subsampled absolute value rectified responses of the first level. In figure 5
we display the first level filters, and the second level filters up to shifts
of size 8 and sign changes of the first level filters.

Figure 5: First and second layer filters from natural images

In order to show the dependence of the filters on the number of filters used,
in figure 6 we display an $8$, $16$, and $64$ element $16\times 16$ dictionary
trained on the same set as above.

Figure 6: Dictionaries with varying numbers of elements trained on natural
images.

## References

  * [1] M. Aharon, M. Elad, and A. Bruckstein.  K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation.  IEEE Transactions on Signal Processing, 54(11):4311–4322, 2006\.
  * [2] Francis R. Bach.  Consistency of the group lasso and multiple kernel learning.  Journal of Machine Learning Research, 9:1179–1225, 2008.
  * [3] Sven Behnke.  Discovering hierarchical speech features using convolutional non-negative matrix factorization.  In IJCNN, pages 7–12, 2008.
  * [4] Alfred M. Bruckstein, David L. Donoho, and Michael Elad.  From sparse solutions of systems of equations to sparse modeling of signals and images.  SIAM Review, 51(1):34–81, 2009.
  * [5] Y. Boureau K. Gregor M. Mathieu Y. LeCun K. kavukcuoglu, P. Sermanet.  Learning convolutional feature hierarchies for visual recognition.  Advances in NIPS, 2010.
  * [6] Stephane Mallat and Zhifeng Zhang.  Matching pursuit with time-frequency dictionaries.  IEEE Transactions on Signal Processing, 41:3397–3415, 1993.
  * [7] Graham Taylor Matthew Zeiler, Dilip Krishnan and Rob Fergus.  Hierarchical convolutional sparse image decomposition.  In The Twenty-Third IEEE Conference on Computer Vision and Pattern Recognition, San Francisco, CA, June 2010.
  * [8] B. Olshausen and D. Field.  Sparse coding with an overcomplete basis set: A strategy employed by v1?, 1997.
  * [9] Ming Yuan, Ming Yuan, Yi Lin, and Yi Lin.  Model selection and estimation in regression with grouped variables.  Journal of the Royal Statistical Society, Series B, 68:49–67, 2006\.