Is Betelgeuse Really Rotating? Synthetic ALMA Observations of Large-scale Convection in 3D Simulations of Red Supergiants

To assess whether a nonrotating RSG can show a dipolar radial velocity map in the ALMA band, we need global 3D RSG models that simulate the multiscale 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π sphere (Goldberg et al. 2022). The CO5BOLD RSG simulations have been extensively used to interpret spectrophotometric, 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) such that it accounts for the heating and cooling effect of the radiation flux. 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_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_pot. For the gravity experienced by the simulation, the stellar mass is M_pot because the self-gravity of the envelope is neglected. However, for actual stars, the total stellar mass would be M_pot + M_env. For detailed descriptions of the general setup, see Chiavassa et al. (2011a), Freytag et al. (2012), and Chiavassa et al. (2024).

We use the 3D radiation hydrodynamic simulations of RSG envelopes presented in Ahmad et al. (2023, and references therein). We chose their model st35gm04n37 for the discussion presented in the main text of this paper (hereafter referred to as our fiducial model or model A). The assumed surface gravity and effective temperature in this model are close to the values observed for Betelgeuse (see Table 2 in Appendix D), but the model is not a perfect match. Therefore, we also analyze an alternative model, st35gm03n020 from Ahmad et al. (2023), hereafter model B, which is more massive. The mass adopted in this model is closer to the inferred mass of Betelgeuse, but it may nevertheless be less appropriate than model A. Arroyo-Torres et al. (2015) show that model B has less extended atmosphere than observed. This is probably because the radiation pressure is neglected in the simulations (Chiavassa 2022). The lower mass and lower surface gravity in model A partly compensate for this. Recent studies therefore used model A to compare with observations of other RSGs that are similar to Betelgeuse (Kravchenko et al. 2019; Climent et al. 2020; Chiavassa et al. 2022). We only present model A in the main text for conciseness but refer the reader to Appendix C.4 for a discussion of the limitations and Appendix D for an analysis of model B. The quantitative results differ, but for both models, we find the same main conclusion that convective motions can mimic rapid rotation at the km s⁻¹ level.

We postprocess the 3D simulations to create synthetic observations for ALMA. The steps toward the synthetic images are as follows.

• 1.  
  Calculate the abundances of atoms and molecules in each cell of 3D simulations to get emissivity and opacity.
• 2.  
  Solve the radiative transfer equations to obtain images of intensity maps.
• 3.  
  Convolve the resulting images with the ALMA beam.

Here we describe the main assumptions. Details and tests for each step are provided in Appendix C. The postprocessing 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 FASTCHEM2 ⁷ (Stock et al. 2018, 2022). FASTCHEM2 has been widely adopted and tested against observations of exoplanetary atmospheres, in particular for ultrahot Jupiters, where the dayside conditions are similar to cool stars (e.g., Kitzmann et al. 2018, 2023). Here, we present the first application of FASTCHEM2 to stellar atmospheres. Chemical equilibrium may not be a good approximation for the shock-dominated atmosphere considered here (e.g., Cherchneff 2006). However, since this study mainly focuses on the velocity field rather than the abundances, the impacts on the main results are expected to be limited.

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 ²⁸SiO lines (vibrational level ν = 2, rotational transition J = 8−7 as observed by ALMA). Detailed equations and opacity sources can be found in Appendix C.2. Since one actual ALMA image of Betelgeuse only needs less than 1 hr to take (K18), the surface motions are approximately frozen during the exposure time of one image. Therefore, no extra integration in time is needed to produce synthetic maps from simulation snapshots.

We assume that the simulated star is located at a distance away from the observers such that the radio photospheric radius approximately matches the ALMA observations in K18 (for the associated limitations, see Appendix C.4). We interpolate the intensities onto a 101² 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{radio}}=\sqrt{\int I(r,\theta ){rdrd}\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, θ) on the 2D image, measured in the spectral window centered at 343.38 GHz as done in the ALMA observations. Its value at the photocenter is denoted as I_star.

To compare our simulations with the observed radial velocity maps, we apply similar analysis methods as in K18 to our generated synthetic observations. 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_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_sys-subtracted radial velocity map within R_radio (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 synthetic radial velocity maps derived from emission lines and not absorption lines, in accordance with the actual analysis of K18.

In Figure 2, we show a selected example of a synthetic ALMA image from our nonrotating RSG simulation and compare it with the ALMA observations of Betelgeuse in K18. The left-hand panels ((a) and (d)) show the original simulation. The middle panels ((b) and (e)) show our synthetic observations after convolution with the ALMA beam. The right-hand panels ((c) and (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 2(a)). The convective motions of ±30 km s⁻¹ can be seen in the radial velocity map of the original simulations (Figure 2(d)). 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 2. For the intensity map, a single hot spot emerges (Figure 2(b)), similar to what is observed by ALMA (Figure 2(c)). 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). Quantitatively, the intensity in the mock image is half of the observed value. This corresponds to a discrepancy of 15% in temperature since the intensity is proportional to temperature to the fourth power. This is mostly due to the fact that the surface temperature is lower and the atmosphere is less extended than observed. 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 km s⁻¹ (Figure 2(d)) into a map with radial velocity variations up to 8 km s⁻¹ (Figure 2(e)) 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 (Figures 2(e) 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 ± 0.10 km s⁻¹ and a residual velocity dispersion of 1.44 km s⁻¹. Following the same procedure for our synthetic image, we obtain very similar values; see Table 1. The resulting ${\chi }_{\mathrm{red}}^{2}$ is larger than observed but comparable within an order of magnitude. As shown in Figure 5 later in the text, if we have underestimated the overall smearing effect of the observational pipeline, the ${\chi }_{\mathrm{red}}^{2}$ may be vastly decreased. In principle, this can be done in the future with the CASA package ⁸ to more properly take into account the antenna configuration, noises, etc.

Table 1. Comparison of Fitting Parameters of a Nonrotating RSG Simulation with ALMA Observations of Betelgeuse

Inferred Parameters	$v\sin i$	${\sigma }_{\mathrm{res}}$	${\chi }_{\mathrm{red}}^{2}$
	(km s−1)	(km s−1)
Mock observation	5.51 ± 0.21	1.68	94.6
(this work, nonrotating)
ALMA observation	5.47 ± 0.10	1.44	20.7
(K18)

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

Download table as:  ASCIITypeset image

Note, however, that our simulation is based on a nonrotating star. The dipolar velocity field shown in the synthetic image is not the result of rotation but of the alignment of the group motions of the surface convective cells. The example shown here illustrates how turbulent motions in the atmosphere of a nonrotating 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.

In the previous section, we showed an example of how a nonrotating simulation can appear to be rotating as a result of underresolved 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 six faces and 80 snapshots (uniformly taken across the 5 yr time span of the relaxed simulation) and plot the distribution in Figure 3.

Zoom In Zoom Out Reset image size

Figure 3 shows that the apparent $v\sin i$ distribution inferred for our nonrotating models peaks at ∼5 km s⁻¹. The observed values for Betelgeuse (gray 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 nonrotating RSG to be interpreted as rotating faster than ∼5 km s⁻¹, the probability still remains as high as ∼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 ∼90% single-epoch observations would show signs of rotation faster than ∼2 km s⁻¹. By visual inspection, the majority of those with inferred $v\sin i\lesssim 2\,\mathrm{km}\,\ {{\rm{s}}}^{-1}$ do not show an obvious dipolar structure that could be mistaken for rotation. When analyzing an alternative simulation with slightly different stellar parameters (see Appendix D, in particular Figure 10), we also find a very high chance of inferring rotation at a few km s⁻¹.

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

Support for the rotation hypothesis (1) 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 yr later.

This interpretation of the HST data is not without controversy. The radial velocities vary between observations taken a few months apart (see Figure 6 in Harper & Brown 2006). Other HST data analyzed by Lobel & Dupree (2001) instead show evidence for a reversal of velocities, and Jadlovský et al. (2023) found no sign of rotation. These findings are more consistent with hypotheses (2) and (3), 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 Section 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). However, so far Betelgeuse is the only RSG with a marginally resolved velocity map available from ALMA.

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

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 1), 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 2) 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 3) 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 may not display a dipolar feature. The error bars quoted by K18 for the projected rotational velocity (0.1 km s⁻¹, corresponding to a relative error of 2%) and the orientation angle of the rotation axis (3.5°, i.e., 1% of a full circle) are so small that the probability of inferring values within these error bars by chance in the future observations is negligible.

Observations of the radial velocity map at increased spatial resolution will also help distinguish the hypotheses we have outlined here. An increase by a factor of 2 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 3 times better. When ALMA upgrades its current 16 km baseline to 32 km, we can expect an (additional) improvement of a factor of 2 (Carpenter et al. 2020).

In Figure 5, we present mock observations with different spatial resolutions. We expect that an increase of a factor of 2 in spatial resolution (panels (c) and (g)) will already marginally resolve the convective structure in both the radial velocity map and the intensity map and be able to distinguish between all three scenarios. On the contrary, decreasing the resolution by a factor of 2 ((a) and (e)) smears out any asymmetrical feature in the intensity map, and the radial velocity map is completely dominated by a dipolar pattern. As the spatial resolution of the interferometer increases (left to right), the measured radio photospheric radius monotonically decreases from ∼30 to ∼25 mas, while the peak magnitude of radial velocity monotonically increases from ∼6 to ∼30 km s⁻¹. As a result of increasing radial velocity magnitude, the fitted $v\sin i$, ${\chi }_{\mathrm{red}}^{2}$, and ${\sigma }_{\mathrm{res}}$ all increase with higher resolution, among which ${\chi }_{\mathrm{red}}^{2}$ increases by 3 orders of magnitude. If future observations can be done with higher resolution, these are the clear signatures to look for to support our hypothesis, although the differences in the frequency range need to be considered.

Zoom In Zoom Out Reset image size

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 precalculate the number density tables that cover the parameter space of $2.5\leqslant {\mathrm{log}}_{10}(T/[{\rm{K}}])\leqslant 5.5,-30\leqslant {\mathrm{log}}_{10}({ \mathcal R })\leqslant -9$, where ${ \mathcal R }={(P/[\mathrm{dyn}\,\ {\mathrm{cm}}^{-2}])/(T/[{\rm{K}}])}^{4}$, using FastChem2. 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 precalculated FastChem2 table are suitable for RSG simulations, we compare the interpolated values to the values given by a MARCS RSG model in Figure 6. The MARCS code was developed to create 1D atmospheric models for evolved stars (Gustafsson et al. 2008), which have 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.

Zoom In Zoom Out Reset image size

In Figure 7, we show an example of the SiO number density distribution in our simulations by taking a 2D slice along the middle cross section. The 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 infrared photosphere due to recombination. The presence of shocks outside of the photosphere can be seen from the SiO number density slice.

Zoom In Zoom Out Reset image size

In our postprocessing, 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

$$\begin{eqnarray}&&{\tau }_{\nu }(\hat{{\boldsymbol{n}}},{\boldsymbol{r}})\equiv {\int }_{0}^{x(\hat{{\boldsymbol{n}}},{\boldsymbol{r}})}{\chi }_{\nu }{{dx}}^{{\prime} }.\end{eqnarray} \tag{ C1 }$$

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

$$\begin{eqnarray}&&{I}_{\nu }(x)={I}_{\nu }(0){e}^{-{\tau }_{\nu }}+{\int }_{{e}^{-{\tau }_{\nu }}}^{1}{S}_{\nu }({\tau }_{\nu }-{\tau }_{\nu }^{{\prime} })d({e}^{-{\tau }_{\nu }^{{\prime} }}),\end{eqnarray} \tag{ C2 }$$

where S_ν is the source function defined by S_ν ≡ η_ν /χ_ν . I_ν , η_ν , and χ_ν are the intensity, emissivity, and opacity at frequency ν. The intensity I_ν (0) at the origin is taken to be 0.

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_ν is the Planck function. The continuum opacity ${\chi }_{\nu }^{\mathrm{con}}$, dominated by 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 ²⁸SiO (ν = 2, J = 8−7) lines are obtained from the ExoMol database ¹⁰ (Tennyson & Yurchenko 2012; Tennyson et al. 2016) using the results of Yurchenko et al. (2022) and further computed assuming local thermal equilibrium (LTE). For this transition, we did not find available information for the collisional rates; therefore, we refrained from performing detailed non-LTE calculations. The line emissivity and opacity contributed by the transition from level i to level j are computed following standard assumptions:

$$\begin{eqnarray}&&{\eta }_{\nu }^{\mathrm{line},{ij}}=\displaystyle \frac{h\nu }{4\pi }{n}_{i}{A}_{{ij}}{\phi }_{\nu }^{{ij}},\end{eqnarray} \tag{ C3 }$$

$$\begin{eqnarray}&&{\chi }_{\nu }^{\mathrm{line},{ij}}=\displaystyle \frac{h\nu }{4\pi }({n}_{j}{B}_{{ji}}-{n}_{i}{B}_{{ij}}){\phi }_{\nu }^{{ij}},\end{eqnarray} \tag{ 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,

$$\begin{eqnarray}&&{\phi }_{\nu }^{{ij}}=\displaystyle \frac{1}{\delta {\nu }_{{ij}}\sqrt{\pi }}\exp \left[-{\left(\displaystyle \frac{\nu -{\tilde{\nu }}_{{ij}}}{\delta {\nu }_{{ij}}}\right)}^{2}\right],\end{eqnarray} \tag{ C5 }$$

$$\begin{eqnarray}&&\delta {\nu }_{{ij}}=\displaystyle \frac{| {v}_{\mathrm{therm}}| }{c}{\tilde{\nu }}_{{ij}},\end{eqnarray} \tag{ C6 }$$

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

$$\begin{eqnarray}&&\begin{array}{l}{\tilde{\nu }}_{{ij}}={\nu }_{{ij}}\Space{0ex}{0.25ex}{0ex}(1-\displaystyle \frac{\hat{{\boldsymbol{n}}}\cdot {\boldsymbol{v}}}{c}\Space{0ex}{0.25ex}{0ex}),\\ {v}_{\mathrm{therm}}=\sqrt{\displaystyle \frac{2{k}_{B}T}{m}}.\end{array}\end{eqnarray} \tag{ C7 }$$

Here, ν_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 toward the observer is blueshifted and has negative values.

We have tested our algorithm against the line transfer code MAGRITTE (De Ceuster et al. 2020a, 2020b, 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 that 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.

An example of our Gaussian line fit procedure is illustrated in Figure 8. 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).

Zoom In Zoom Out Reset image size

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_sys-subtracted radial velocity map has a velocity distribution of $v(x,y)=v\sin i\times (x\cos \theta +y\sin \theta )/{R}_{\mathrm{radio}}$, where the fitted parameters are the position angle θ and the inferred $v\sin i$, and x, y are the coordinates of the 2D radial velocity map.

We caution that there are different radii involved in interpreting the observations and comparing with 1D or 3D models: the infrared photosphere R_infrared approximately used as outer boundaries in 1D models, the radio photosphere R_radio obtained from ALMA continuum intensity maps, and the molecular shell radius R_shell probed through molecular lines in ALMA. It generally follows that R_infrared < R_radio < R_shell, and for Betelgeuse, each two nearby radii are different by a factor of ∼1.3 (K18). Instead of using the radio photospheric radius R_radio for the fit, we suggest it is more consistent with the physical scenario to use the radius R_shell ≈ 1.3R_radio 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}\,\ {{\rm{s}}}^{-1}$ at the molecular shell. Assuming that the molecular shell corotates with the radio photosphere, this translates to a rotational velocity at the radio photosphere of 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}\,\ {{\rm{s}}}^{-1}$ at the radio photosphere. However, if the molecular shell and the radio photosphere are not fully coupled, e.g., assuming a 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}\,\ {{\rm{s}}}^{-1}$ at the radio photosphere. Since the infrared photosphere R_infrared is yet another factor of 1.3 smaller than the radio photosphere R_radio, the measurement by K18 would give $v\sin {i}_{\mathrm{infrared}}\approx 4\,\mathrm{km}\,\ {{\rm{s}}}^{-1}$ at the infrared photosphere R_infrared assuming constant angular frequency (corotation) outside the star, or $v\sin {i}_{\mathrm{infrared}}\approx 12\,\mathrm{km}\,\ {{\rm{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_radio for the rotation fit, with the underlying assumption that the molecular shell corotates with the radio photosphere.

We place the simulated star at a distance away from the observers such that it has an inferred radio photospheric radius R_radio of ∼30 mas. However, this corresponds to a distance of 130 pc, consistent with the original Hipparcos parallax (${131}_{-23}^{+35}$ pc; ESA 1997; Perryman et al. 1997) and its revised value (${152}_{-17}^{+22}$ pc; Van Leeuwen 2007) but much closer than the distance derived from Hipparcos data combined with radio observations (${197}_{-35}^{+55}$ pc; Harper et al. 2008) and its revised value (${222}_{-34}^{+48}$ pc; Harper et al. 2017). As a comparison, asteroseismic constraints suggested a medium value (${168}_{-15}^{+27}$ pc; Joyce et al. 2020). This discrepancy between the value adopted in the synthetic images and the observed value is partly because the fiducial 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 (Chiavassa 2022).

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 8 compared to Figures 1 and 9 in K18). This could result from the small atmospheric extension of the simulations or not including nonequilibrium chemistry, non-LTE populations, dust emission, maser, or scattering in the postprocessing. 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 height dependence in the velocity field from the surface to the molecular layer. From the physical point of view discussed in Appendix B, the convective velocity v_c scales very weakly to the surface energy flux F and surface density ρ_s as ${(F/{\rho }_{s})}^{1/3}$, so the magnitude of the convective velocity is only weakly affected by the limitations of simulations.

A more uncertain aspect is whether the convective structure in our simulations can represent actual RSG convection. The convective length scales are subject to the uncertainties of the 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). For typical RSGs, where the bottom of the convective envelope is less than 0.05 of the stellar radius, we expect that similar dipolar flows may also dominate in the convective envelope. These dipolar flows are partly 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.