# The Great Dimming of Betelgeuse: The photosphere as revealed by tomography over the past 15 yr*

Daniel Jadlovský<sup>1,2</sup>, Thomas Granzer<sup>3</sup>, Michael Weber<sup>3</sup>, Kateryna Kravchenko<sup>4</sup>, Jiří Krtička<sup>1</sup>, Andrea K. Dupree<sup>5</sup>, Andrea Chiavassa<sup>6</sup>, Klaus G. Strassmeier<sup>3</sup>, and Katja Poppenhäger<sup>3</sup>

Received 5 December 2023 / Accepted 14 February 2024

## ABSTRACT

**Context.** Betelgeuse, a red supergiant star of semi-regular variability, reached a historical minimum brightness in February 2020, known as the Great Dimming. Even though the brightness has returned to the values prior to the Great Dimming now, it continues to exhibit highly unusual behavior.

**Aims.** Understanding the long-term atmospheric motions of Betelgeuse and its variability could be a clue to the nature of the Great Dimming and the mass-loss process in red supergiants. Our goal is to study long-term dynamics of the photosphere, including during the Great Dimming.

**Methods.** We applied the tomographic method, which allows different layers in the stellar atmosphere to be probed in order to reconstruct depth-dependent velocity fields. The method is based on the construction of spectral masks by grouping spectral lines from specific optical depths. These masks are cross-correlated with the observed spectra to recover the velocity field inside each atmospheric layer.

**Results.** We obtained about 2800 spectra over the past 15 yr that were observed with the STELLA robotic telescope in Tenerife. We analyzed the variability of five different layers of Betelgeuse's photosphere. We found phase shift between the layers, as well as between the variability of velocity and photometry. The time variations of the widths of the cross-correlation function reveal propagation of two shockwaves during the Great Dimming. For about 2 yr after the dimming, the timescale of variability was different between the inner and outer photospheric layers. By 2022, all the layers seemingly started to follow a similar behavior as before the dimming, but pulsating with higher frequency corresponding with the first overtone.

**Conclusions.** The combination of the extensive high-resolution spectroscopic data set with the tomographic method revealed the variable velocity fields in the photosphere of Betelgeuse, for the first time in such detail. We were also able to find new insights related to the Great Dimming event and its aftermath, namely the discovery of another shockwave and the subsequent rearrangement of the photosphere. Our results demonstrate that powerful shocks are the triggering mechanism for episodic mass-loss events, which may be the missing component to explain the mass-loss process in red supergiants.

**Key words.** techniques: photometric – techniques: spectroscopic – stars: late-type – stars: mass-loss – stars: winds, outflows

## 1. Introduction

Betelgeuse (α Orionis) is a nearby red supergiant (RSG) star. Its semi-regular variability is characterized by two main timescales common to many RSGs. These are (i) a shorter period of roughly 400 days, attributed to the radial pulsations of the atmosphere in the fundamental mode (FM, Kiss et al. 2006; Chatys et al. 2019; Joyce et al. 2020; Jadlovský et al. 2023); and (ii) a longer period of about 2100 days (5.6 yr), which is of a less certain origin. The longer period has been attributed to the flow timescales of giant convective cells (Stothers 2010), which

López Ariste et al. (2022) found some evidence for through spectropolarimetry imaging. It has also been attributed to processes that drive a similar mode of variability in asymptotic giant branch (AGB) stars, where it is known as the long secondary period (LSP), most likely caused by nonradial, low-degree oscillations due to gravity modes. However, other mechanisms were also proposed: for example, magnetic activity (Wood 2000; Wood et al. 2004). In many RSGs, Kiss et al. (2006) and Chatys et al. (2019) also found higher overtones of the FM period, primarily the first overtone. For Betelgeuse, Joyce et al. (2020) detected a photometric period of about 185 days, which is likely the first overtone mode. Granzer et al. (2022) found a period of about 216 days based on radial velocities as well as evidence for the second overtone. However, for all the modes of variability, the observed length and amplitude of each cycle is not constant (Kiss et al. 2006).

<footnote n="1">¹ Department of Theoretical Physics and Astrophysics, Faculty of Science, Masaryk University, Kotlářská 2, Brno, 611 37, Czech Republic
e-mail: jadlovsky@mail.muni.cz</footnote>

<footnote n="2">² European Southern Observatory (ESO), Karl-Schwarzschild Str. 2, 85748, Garching bei München, Germany</footnote>

<footnote n="3">³ Leibniz-Institut für Astrophysik Potsdam (AIP), An der Sternwarte 16, 14482 Potsdam, Germany</footnote>

<footnote n="4">⁴ Max Planck Institute for Extraterrestrial Physics (MPE), Giessenbachstrasse 1, 85748, Garching bei München, Germany</footnote>

<footnote n="5">⁵ Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, MS-15, Cambridge, MA 02138, USA</footnote>

<footnote n="6">⁶ Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Lagrange, CS 34229, 06304 Nice Cedex 4, France</footnote>

<footnote>* The reduced STELLA echella spectra and determined radial velocity are only available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/685/A124</footnote>
Betelgeuse underwent an unprecedented historical minimum of brightness, the so-called Great Dimming, that culminated in February 2020 when its brightness dropped by $V \sim 1.6$ mag (Guinan et al. 2020). Interferometric observations by Montargès et al. (2021) revealed that the star dimmed asymmetrically, as the southern hemisphere of Betelgeuse was much darker than its northern counterpart. Many theories have been put forward to explain the Great Dimming, such as a significant decrease in temperature or the formation of star spots (Dharmawardena et al. 2020; Harper et al. 2020), an increased molecular opacity in the outer layers (Kravchenko et al. 2021), or a surface mass-loss event (Dupree et al. 2020, 2022). Interferometric observations appear to be compatible with both leading theories: the cold spot on the surface as well as a dust clump above the photosphere (Montargès et al. 2021; Cannon et al. 2023; Drevon et al. 2024).

In the detailed study that combined the previous findings, Dupree et al. (2022) suggested that the Great Dimming was caused by a shockwave (reported by Kravchenko et al. 2021) that formed in the photosphere and resulted in substantial surface mass ejection of matter along our line of sight in September-November, 2019. The ejected material reached the outer atmosphere in about 6 months, eventually creating the molecules and dust as it reached cooler regions. Dupree et al. (2022) and Jadlovský et al. (2023) also observed that following the Great Dimming, the star pulsated with shorter periods than 400 days.

New models by Davies & Plez (2021), which successfully incorporated stellar wind and, for the first time, reproduced extended atmospheres and related observed spectral features of RSGs, generally support the scenario proposed by Dupree et al. (2022), although the authors show that the dimming can be explained by TiO molecule absorption, without the need for dust condensation in the ejected material. Nonetheless, Taniguchi et al. (2022) reported an increased extinction by oxygen-rich dust that contributed to the Great Dimming.

In the present paper, we intend to reveal long-term atmospheric motions of Betelgeuse and its variability by applying the tomographic method by Alvarez et al. (2001a) and Kravchenko et al. (2018) to a time series of high-resolution spectra of Betelgeuse. In Sect. 2, we describe our data set and the application of the tomographic method. In Sect. 3, we present the results. In Sect. 4, we take a more detailed look at the Great Dimming and its aftermath, based on our data. In Sect. 5, we discuss the shocks and excitation of higher overtones, and we give our conclusions in Sect. 6.

# 2. Observations and methodology

## 2.1. Observations

We employed spectra observed with the Stellar Activity (STELLA) echelle spectrograph (SES) mounted on two fully robotic 1.2 m telescopes at the Izanã Observatory in Tenerife, which is operated by AIP (Strassmeier et al. 2004, 2010). The resolving power is $R \approx 55\,000$ with a three-pixel sampling per resolution element. The spectra typically have signal-to-noise ratios (S/N) of 100–600. Exposure times range from five to 20 seconds. The spectra cover the 390–880 nm wavelength range. As of January, 2024, the full time series of Betelgeuse consists of about 2800 spectra (about 6000 exposures in total) obtained during the years 2008–2024. Most of these observations were originally presented by Granzer et al. (2021, 2022), where the telescope, observations, and methodology are also briefly described.

## 2.2. Methodology

### 2.2.1. Envelope tomography

The tomographic method was first introduced by Alvarez et al. (2000, 2001a). The method allows us to probe different layers in the stellar atmosphere and recover the corresponding velocity fields. This is accomplished by sorting spectral lines into groups (masks) based on their formation depth. The masks are cross-correlated with observed or synthetic stellar spectra and provide velocity fields within different atmospheric layers.

Cross-correlation (CCF) functions of some Mira-type stars appear double-peaked during the maximum brightness (Alvarez et al. 2000). One of the main goals of using the tomographic method was to show that this feature can be explained by the passing of shockwaves through the atmosphere, which is called the Schwarzschild scenario (Schwarzschild 1952); that is, the second peak would correspond to the rising shock front. Alvarez et al. (2000, 2001a,b) successfully applied the method and recovered the Schwarzschild scenario on spatial and temporal scales in several Mira-type stars. The only limitation of the method was that it assumed that the line depression originates in layers where the optical depth $\tau_{\lambda}$ is $2/3$. However, that was shown not to be the case for weak lines (Albrow & Cottrell 1996). It was also difficult to make quantitative predictions as the geometric radius associated with each mask was not known.

Kravchenko et al. (2018) improved the method by including the computation of the contribution function to calculate the formation depth of spectral lines. The method was then validated on a 3D radiative-hydrodynamics simulation of an RSG atmosphere (Chiavassa et al. 2011). Kravchenko et al. (2019) applied the tomographic method to high-resolution spectroscopic time-series observations of the RSG star $\mu$ Cep to interpret its photometric variability. Finally, Kravchenko et al. (2020) validated the tomographic method on spectro-interferometric VLTI/AMBER observations by recovering a link between optical and geometrical depth scales based on comparison to dynamical model atmospheres.

Afterward, Kravchenko et al. (2021) also applied the tomographic method to high-resolution spectroscopic time-series observations of Betelgeuse taken with the HERMES spectrograph from 2015–2020 to interpret its photometric variability, including the Great Dimming event. The resolution of HERMES is 86 000, which is larger than resolution of STELLA, and their data set included about 30 spectra. They found two shockwaves preceding the dimming (in February, 2018 and January, 2019), the former amplifying the latter. They proposed that the succession of these two shocks resulted in an outflow (as shown by velocity gradients) and an increase of molecular opacities as their interpretation of the dimming. They did not rule out the dust scenario, but they could not confirm it based on spectroscopic observations.

To apply the tomographic method in this paper, we used a spectral template<fnref n="1" /> (spherical symmetry, $T_{\text{eff}} = 3400$ K, $\log g = -0.4$, $v_{\text{microturb}} = 2$ km s⁻<fnref n="1" /> and solar metallicity) calculated with the Turbospectrum synthesis code (Plez 2012) with MARCS models (Gustafsson et al. 2008) and VALD3 atomic and molecular line lists (Ryabchikova et al. 2015). Furthermore, to create tomographic masks that correspond to different layers of Betelgeuse’s photosphere, we used the same masks as Kravchenko et al. (2019). Therefore, the spectral template for Betelgeuse was modified to include only the lines from each tomographic mask; that is, we created five new templates with

<footnote n="1">¹ Downloaded from https://marcs.oreme.org</footnote>
<table><caption>**Table 1.** Lists of spectral lines from Kravchenko et al. (2019) used in each mask, $\tau_0$ limits, and numbers of spectra, which were left from the original data set after successful cross-correlations and subsequent cleaning.</caption><thead><tr><th>Mask</th><th>log τ<sub>0</sub> limits</th><th>No. of lines</th><th>No. of spectra</th></tr></thead><tbody><tr><td>C1</td><td>-1.0 &lt; log τ<sub>0</sub> ≤ 0.0</td><td>135</td><td>1679</td></tr><tr><td>C2</td><td>-2.0 &lt; log τ<sub>0</sub> ≤ -1.0</td><td>1750</td><td>2290</td></tr><tr><td>C3</td><td>-3.0 &lt; log τ<sub>0</sub> ≤ -2.0</td><td>1199</td><td>2309</td></tr><tr><td>C4</td><td>-4.0 &lt; log τ<sub>0</sub> ≤ -3.0</td><td>433</td><td>2311</td></tr><tr><td>C5</td><td>-4.6 &lt; log τ<sub>0</sub> ≤ -4.0</td><td>378</td><td>2312</td></tr></tbody></table>

selected lines. This was achieved by removing all lines from the template that were not included in the masks (for each mask). This is done differently than in Kravchenko et al. (2019, 2021), where Dirac distribution was used as line profiles.

Table 1 lists the masks used in this paper. We used the same masks as Kravchenko et al. (2021), except for mask C1, where we had to use fewer lines in order to increase the signal. Lines in mask C1 have the highest optical depths, while each following mask has lower values of optical depths than the one before it, all relative to a reference optical depth $\tau_0$ computed at 5000 Å, as described in Kravchenko et al. (2018). Therefore, mask C1 corresponds to the innermost layer of the probed photosphere and mask C5 corresponds to the outermost layer.

### 2.2.2. Cross-correlation functions

The data were reduced using the IRAF-based pipeline SESDR in its version 4.0 (Weber et al. 2008, 2012). All spectra were corrected for echelle blaze and wavelength calibrated with con-secutively recorded Th-Ar spectra. To derive cross-correlation functions (CCFs), the radial velocity determination is done order-by-order, that is, each echelle order is cross-correlated with the mask separately. Afterward, the cross-correlation func-tions are summed up, and the radial velocity is determined by a Gaussian fit (for a single observation).

Many of the spectra were removed due to failed guiding of the telescope or other various instrument-related issues (about 7%). The numbers of successfully cross-correlated spectra left after the cleaning are listed in Table 1.

## 3. Results

This section is focused on characterizing the overall variabil-ity of Betelgeuse based on the application of the tomographic method to our data set. Our results are also compared to those of Kravchenko et al. (2021), who applied the tomographic method to Betelgeuse for the first time.

### 3.1. Brightness

The first panels of Figs. 1 and 2 show the evolution of brightness adopted from the American Association of Variable Star Observers (AAVSO), Solar Mass Ejection Imager (SMEI, Jackson et al. 2004; Hick et al. 2007), and Ogane Hikari Obser-vatory (Ogane et al. 2022). The first two data sets consist of more data than plotted in the figures. For SMEI, a cleaning algorithm as described in Paunzen et al. (2021) and Jadlovský et al. (2023) was applied.

The combination of these data sets covers the semiregular photometric variability of the star and the Great Dimming well. Both main modes of variability are clearly visible in the data, as well as variability of a shorter timescale than the FM period.

### 3.2. Relative temperature changes

The line-depth ratio of V I (6251.83 Å) and Fe I (6252.57 Å) is sensitive to temperature changes due to excitation potentials of corresponding levels; therefore, it can be used as a proxy for tem-perature. The ratio becomes larger for decreasing temperature in the photosphere (Gray 2000, 2008). Calibrating these ratios to absolute values of effective temperature would be difficult. How-ever, the relative variability of temperature is sufficient for us, as we primarily use it for plotting hysteresis loops in Sect. 3.6.

The ratios are plotted alongside the brightness in the upper panel of Figs. 1 and 2. The derived relative temperature changes are mostly in phase with the photometric variability, similarly to the derived effective temperature in previous studies (e.g., Kravchenko et al. 2019; Dupree et al. 2022).

### 3.3. Radial velocities

The second panels of Figs. 1 and 2 reveal kinematics of five different layers of the photosphere using five cross-correlation tomographic masks from Table 1. Thanks to the abundance of observations, we are able to see the long-term evolution of dif-ferent layers of Betelgeuse in such a detail for the first time. The zero-point is set at 20.7 km s⁻¹; this corresponds to the center-of-mass (CoM) velocity of Betelgeuse (the average value between Harper et al. 2017 and Kervella et al. 2018).

The usual maximal peak-to-peak amplitudes of velocity vari-ations in all photospheric layers reach up to ~5 km s⁻¹, but often the amplitudes are smaller, especially when there is a variability within shorter timescales than the FM period. All these modes of variability are also modulated by the LSP period, which is clearly evident in the data and is represented in Figs. 1 and 2 by a sinusoidal curve. This curve was determined based on a Lomb-Scargle periodogram performed on the radial velocity data (Granzer et al. 2022). Parameters of the sinusoidal curve are a period of 2081 days, an amplitude of 2.54 km s⁻¹, and initial phase of 0.71. The shorter periods of variability, likely the higher overtones of the FM period, appear to become more prominent when the photospheric velocities are lower than the velocity of CoM of the star for several cycles of the FM period. This is shown in Fig. 1 from 2008–2010 and in Fig. 2 from 2020–2023, and it is further discussed in Sects. 4 and 5.

When rising material starts to propagate through the photo-sphere, the innermost layer, C1, is the first to be affected, and it is mostly in phase with brightness V and temperature variabil-ity; namely, the minima occur in approximately the same point of time. The time shift of velocity minima between the inner-most (and thus also the brightness and temperature minima) and outermost layers can be up to about 50 days (10 days between neighboring layers), which is quite common from 2017–2023, and the overall scenario is well illustrated by the minimum at the beginning of 2019. The values of time shifts correspond to the ones reported by Kravchenko et al. (2021). However, from 2008–2017, the time shifts are much smaller, only up to about 10 days in total, while the velocity shifts between the layers are also less significant. Nonetheless, we have fewer minima covered by the data in this time period. During the velocity maxima, the time and amplitude shifts are also shorter. Because the (observed) length and amplitude of each cycle is not constant, we only
<figure><img src="image_1.png" /><figcaption>**Fig. 1.** Variability of Betelgeuse from 2008–2016 *First panel:* light curve comprising of data from AAVSO, Ogane, and SMEI. Temperature is also included via the line-depth ratio of V I and Fe I. *Second panel:* radial velocity plot determined from STELLA spectra. The spectra were cross-correlated with five masks corresponding to five different layers of photosphere. The same color-coding is applied in the fourth panel. The velocity of the center of mass (20.7 km s⁻¹) was subtracted. Both main periods, the FM and LSP, are clearly evident in the data, while we can also see the presence of higher overtones of the FM period between years 2009 and 2010. *Third panel:* gradients of radial velocity from masks in second panel. At most times the gradients are not very high. *Fourth panel:* FWHMs of CCFs. The FWHMs are correlated with the variability of the star, to a certain extent.</figcaption></figure>

give approximate values of time and amplitude shifts between
individual radial velocity curves.

**3.4. Velocity gradients**

The third panels of Figs. 1 and 2 show curves of velocity gradi-
ents, that is, the evolution of velocities relative to the outermost
layer C5. The gradients were calculated in the same manner as in
Kravchenko et al. (2021); that is, we defined gradients for layers
C1–C4 as vC5 − vCi, where i = 1, 2, 3, 4. Thus, when the gradients

are negative (v<sub>C5</sub> − v<sub>Ci</sub> < 0), they signify the expansion of regions
over which the gradient was taken, whereas when the gradients
are positive (v<sub>C5</sub> − v<sub>Ci</sub> > 0), they signify contraction.

The gradients are quite small at most times, within 2 km s⁻¹
of the zero point, and quite often even at the zero point. In other
words, at most times all the layers are moving with nearly the
same velocity, such as between the years 2014 and 2017.

In cycles when the velocity variability is stronger, the gradi-
ents reveal to us that when the material is rising, the inner layers
of the photosphere are typically being contracted (see, e.g., the
<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Same as Fig. 1, but for 2016–2024. *First panel:* light curve shows historical Great Dimming. *Second panel:* prior to Dimming, we can see a major long-lasting outflow velocity, followed by extreme infall velocity. *Third panel:* during Great Dimming, gradients were uniquely negative, suggesting a strong expansion of the inner layers. *Fourth panel:* before and after Dimming, there are two significant peaks in the FWHMs, while the latter peak was much stronger. Both peaks are caused by a passing shockwave.</figcaption></figure>

positive gradients at the beginning of years 2013 and 2019 in Figs. 1 and 2, respectively). On the other hand, when the material in the layers is infalling, all the layers usually move with a similar velocity. Except during the Great Dimming, we only rarely see a significant expansion of the atmospheric layers.

### 3.5. FWHMs

The bottom panel of Figs. 1 and 2 shows the full width at half maximum (FWHM) of the CCFs. The average FWHM in layer C1 is about 17 km s⁻¹. The FWHM increases in the upper layers, up to about 37 km s⁻¹ in layer C5. Our CCF profiles resemble

the ones by Kravchenko et al. (2021), which reported secondary components shifted by 10–15 km s⁻¹ with respect to the main peak. We have also found many similar weak asymmetries, especially in CCFs of mask 1. A sample of our CCFs showing the typical asymmetries is plotted in Fig. 3.

Similarly to Kravchenko et al. (2021), we find the FWHM of CCF to be a reliable indicator of a passing shockwave. We can see a significant increase (by a factor of about 1.3) in FWHMs at the beginning of 2019 (Fig. 2, also shown in Kravchenko et al. 2021), followed by an even more unprecedented increase (by a factor of about 1.5 in outer layers) at the beginning of 2020. Meanwhile, in the rest of the data, the variability of the FWHMs
<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Preview of CCF profiles sampled during cycles related to Dimming. The values of CCF intensity (y-axis) are inverted. Vertical red lines show radial velocity determined. In CCFs of mask C1, we can often clearly see that the central components are asymmetric, suggesting a secondary component. Other inner layers sometimes also show weak asymmetries.</figcaption></figure>

is considerably smaller. This makes it clear that the Great Dimming was truly an extraordinary and unprecedented event. While the first shockwave in the beginning of 2019 was already reported by Kravchenko et al. (2021), the stronger shockwave in 2020 has not yet been reported in any study.

For both shockwaves, we can see a time shift of the emergence of a shockwave in a similar manner to that seen in the RVs; that is, the shockwave is moving from the innermost mask to the uppermost mask. However, there are some differences between the shockwaves in 2019 and 2020. During the propagation of the first shockwave, the FWHM peak of the innermost layer C1 arises in a similar point of time as its velocity minimum, and the peaks of outer masks are increasingly more time shifted in relation to layer C1. In total, the time shift between layers C1 and C5 is about 75 days. Likewise, we can see a similar time shift during the shockwave in 2020, even though the FWHM peaks are larger in outer layers compared to the previous shockwave, while in the innermost layers (C1 and C2) the peaks are barely visible at all, and in layer C3 the peak is still smaller than during the previous shock. This suggests that the initial low-amplitude shock steepened significantly as it propagated through the photosphere. The shock was amplified by a factor of 1.5 by the time the shock front reached photospheric layers with a reference optical depth $log \tau_0$ smaller than -3.0.

The shockwave in February 2018 reported by Kravchenko et al. (2021) does not seem to be represented in our data by any significant FWHM peak (at least not one comparable to the two main peaks), although there is definitely a small increase of FWHM in 2018. That suggests that this shockwave was of a much smaller amplitude than the two that accompanied the dimming. Supposing this was indeed a shock of a smaller scale, then similar shocks of weaker amplitudes are quite frequent, as in Figs. 1 and 2 we can see many peaks at this variability scale.

Additionally, we also remark that resolution of HERMES is about 86 000, which is much higher than that of STELLA. Therefore, it is expected that our instrument would be less sensitive to asymmetries and possible line-doubling events.

### 3.6. Hysteresis loops

The phase shift between variations of temperature and velocity results in hysteresis loops, that is, the evolution of velocity and temperature during a (photometric) cycle. Gray (2008) was able to construct hysteresis loops based on observations of Betelgeuse and proposed the following interpretation. During a regular pulsation cycle, the material starts heating (for the beginning of a cycle defined in the last temperature minimum), while the velocity is constant. Eventually, once the temperature is high enough, the material starts to rise up in the atmosphere. When the material is rising, it begins to cool, and hence it starts to infall again after some time. However, this interpretation assumes a stationary convection, and it thus cannot explain the observable time-dependent effects: specifically the phase shift between the temperature and velocity variations, as demonstrated by Kravchenko et al. (2019, 2021). Based on simulations and theoretical predictions, Kravchenko et al. (2019) investigated the role of convection and acoustic waves in driving the hysteresis loops. They showed that the timescales of hysteresis loops are similar to the sound-crossing timescales in outer layers. Therefore, they proposed that hysteresis loops have an acoustic origin. All these papers also showed that the loops turn counterclockwise and that timescales of loops are similar to the FM period.

We constructed hysteresis loops for Betelgeuse, and it appears that the behavior described in the aforementioned papers can generally be observed in many of the cycles, even though the exact shapes and amplitudes of each loop are often quite different from each other. The results are available in Appendices A.1–A.3. Cases when a clear determination of the shape or length of a loop was not possible originate in the lack of data (typically in older cycles) and poor data sampling. Additionally, the determination of the onset of cycles (based on photometric minima) is not very precise due to the complex variability of the star.

Notwithstanding all these issues, it is quite apparent from the data that the oscillations shorter than the FM period rarely
<figure><img src="image_4.png" /></figure>

**Fig. 4.** Detailed view of Great Dimming and its aftermath from Fig. 2. The upper panel shows variability of brightness and temperature; the lower panel shows variability of velocity. The first red arrow denotes the unusual outflow that lasted ~200 days. During this outflow, the ejection of material as proposed by Dupree et al. (2022) could easily happen. What followed were extreme infall velocities, while the maximum infall velocities of inner layers were higher than those of upper layers, which is unprecedented in the data. In the aftermath of the dimming and the second shockwave, the outer and inner layers behaved differently. Inner layers (C1–C3) exhibited a separation of local outflow maxima corresponding to the first overtone (the third and fourth red arrows), whereas the outer layers appeared to have a substantially longer separation, likely corresponding to the FM. Since the beginning of 2022, the layers appear more synchronized with each other again, but the dominant mode of pulsations is now the first overtone.

form full hysteresis loops. Usually, such curves seem to cover about half of a loop. On the other hand, the FM cycles seem to form full loops. Therefore, this suggests that the timescales of hysteresis loops are similar to those of the FM period.

Furthermore, in Appendix A.4 we show the evolution of velocity and temperature for the LSP cycles. This is done by binning the velocity and relative temperature changes in 400 d intervals, which approximately corresponds to an average length of the FM period. Quite surprisingly, the material in the photosphere also appears to form hysteresis loops during this timescale, superimposed over the shorter loops, although it is not certain whether the first LSP cycle forms a full loop. These longer loops could correspond to the turnover of material in giant convection cells on the surface.

# 4. Great dimming

## 4.1. Emergence

We revealed several unique features that appeared during the Great Dimming (a closer look in Fig. 4). Firstly, there was the powerful shockwave (Kravchenko et al. 2021) that first appeared in the innermost layer, C1, at the beginning of 2019 and then continued to propagate upward. Simultaneously, it caused a strong outflow of rising material (or enhanced it) that reached its maximum by the middle of 2019, with a peak-to-peak amplitude of about 7 km s⁻¹ in some masks. Although the outflow velocity during the peak was quite similar to some of the preceding

maxima, what made this event special was that the outflow velocity stayed at near maximum values for ~200 days (the first red arrow in Fig. 4). This is when material of higher density would be ejected from the atmosphere as reported by Dupree et al. (2020, 2022). Granzer et al. (2022) reported that the pulsations have already been excited to higher overtones during this time.

What followed was a massive change of velocity; the velocity of layers C1-C4 changed quite abruptly by about 10 km s⁻¹ between October, 2019 and April, 2020, while for layer C5 it was by about 8 km s⁻¹. Most importantly, the inner layers were infalling much faster than the layer C5, causing the inner layers to expand dramatically (as shown by the extreme gradients). When the maximum infall velocity was reached by April, 2020, the inner layers, especially C2-C4, abnormally (compared to other cycles) had greater infall velocity than layer C5.

In summary, the sequence of events that led to the Great Dimming began at the beginning of 2019 with the appearance of the shockwave in the photosphere (Kravchenko et al. 2021), followed by an enhanced chromospheric flux in September-November 2019 (Dupree et al. 2022). The brightness reached its historical minimum in February 2020 (Guinan et al. 2020), followed by a minima of photospheric velocities in April, 2020 (Granzer et al. 2022), and then also by a maximum outflow velocity of stellar wind at the end of March (Jadlovský et al. 2023). Meanwhile, the brightness of Betelgeuse quickly began to restore its normal values. Regardless, the events that followed the dimming also continued to be unprecedented.
## 4.2. Aftermath

By the time the innermost layers reached maximum infall velocity during the Great Dimming, a rather small peak appeared in the FWHMs of C1-C2, possibly as the origin of the next shock-wave, which eventually became even stronger than the first one (based on FWHMs behavior). By the time mask C5 reached its maximum infall velocity, the shockwave already appeared to be passing through mask C3, but the height of the FWHMs peak was still smaller than in the same mask during the previous shockwave. Simultaneously, the layers were starting to contract again. By the time the velocities of all layers were rising again, the shockwave passed through layers C4 and C5 in April-May 2020. The passing shockwave appears to reach its maximum FWHM in layer C5 in the same time the data coverage ends. In these two layers, the shockwave was much stronger than the previous one, and the strongest in the entire data set.

Dupree et al. (2022) reported that the atmosphere was left less dense following the surface mass ejection. That could help explain the appearance of the new shock directly after the Great Dimming, and especially its amplified strength. However, the results by Dupree et al. (2022) show that unlike during the previous event, this shock was not followed by a comparable enhancement of chromospheric flux. There was only a slight enhancement of the flux in August, 2020.

Either way, the resulting outflow seems to have one of the highest peak-to-peak amplitudes in the data. The amplitude of the peak cannot be determined precisely, as the peak happened sometime between May and August, 2020, a time period not covered by our observations. Fortunately, the photometric data by Taniguchi et al. (2022) reveal that the maximum brightness was reached at approximately the moment our data coverage ends, just after the beginning of May, 2020. Thus, the maximum outflow velocity was probably reached by the middle of the observation gap, which suggests that by that time, the rising material was able to reach quite high values of outflow velocity. It seems this outflow velocity maximum followed about 250 days after the previous one, which is shorter than the FM period.

The second shock led to a remarkable divergence of the inner (C1-C3) and outer (C4-C5) layers of the atmosphere, which is shown in detail in Fig. 4. The outer layers continued their previous movement, while the inner layers were oscillating with shorter modes of variability of 100-200 days, likely the first or other low-order overtones. This unprecedented dissonance continued for about a year and half. By 2022, all the layers of the photosphere were finally synchronized, in the same way as before the dimming. However, this time the shorter mode of pulsations of about 200 days was dominant (Dupree et al. 2022; Jadlovský et al. 2023; Granzer et al. 2022), which is very likely the first overtone.

However, at the beginning of 2023, the gradients show that inner layers were being significantly expanded (the second strongest expansion after the dimming). This suggests that another major event might have occurred. In April, 2023, this was followed by another brightening. Furthermore, over the last six months (since July, 2023), the brightness surprisingly shows no significant variations, while the temperature and radial velocity continue its variability.

The LSP of about 2100 days appears to remain active, but whether the FM period regains its dominance remains to be seen. Based on the previous behavior of the star, it seems likely that the FM shall indeed return. Hydrodynamic simulations by MacLeod et al. (2023) suggest that the first overtone will

eventually be reduced by the convective damping in about 5 yr (post-dimming).

# 5. Discussion

## 5.1. Shockwaves

The spectral lines considered here come from higher up in the photosphere where the reference optical depth is smaller than one. In these layers of RSGs, the photospheric activity is characterized by interactions between nonradial pulsations and shockwaves (Chiavassa et al. 2010). Based on 3D radiative-hydrodynamics simulations of RSG stars (Chiavassa et al. 2011, 2024; Kravchenko et al. 2019) and AGB stars (Freytag et al. 2017; Liljegren et al. 2017; Chiavassa et al. 2024), the shockwaves originate from acoustic waves produced in stellar interior by convection. As the wave hits the surface, it is slowed down and compressed due to the drop in temperature and sound speed. It increases its amplitude due to the decrease in density and thus turns into a shock, propagating to the outer atmosphere. The density is lower in outer layers, which causes the shockwaves to become even stronger. Consequently, the influence on the photosphere would also become more significant. As we can see in the bottom panel of Fig. 2, both major shockwaves indeed become stronger as they reach the upper layers.

Mira-type stars follow the Schwarzschild scenario (as described in Sect. 2.2.1); that is, evident line doubling occurs in the CCFs, revealing the upward motion of the shock front, when a shockwave is passing through an atmosphere (Alvarez et al. 2000). However, as shown by Josselin & Plez (2007) and Kravchenko et al. (2019), the CCF profiles of RSGs do not clearly follow the Schwarzschild scenario. Instead, the secondary components primarily cause CCFs to become asymmetric, as in these stars the motions within the atmosphere are more complex. The asymmetries are related to supersonic velocity fields due to the motion of giant convective cells and possibly also non-spherically symmetric shockwaves. Based on 3D radiative-hydrodynamics simulations, Kravchenko et al. (2019) demonstrated that the intensity of secondary components is directly related to the emitting surface on the star. For example, when a significant fraction of the surface is rising, a blue secondary component appears in a CCF. Additionally, Ma et al. (2024) demonstrated that large-scale convective motions can produce components that may be mistaken for the signs of rotation.

Similarly to Kravchenko et al. (2021), we can see in our CCFs that even during the dimming event there is no obvious hint of the Schwarzschild scenario in Betelgeuse, and rather than visible line-doubling events we just observe weak asymmetries (which may be due to resolution). Nonetheless, the FWHMs demonstrate that two powerful shockwaves appeared in Betelgeuse before and after the dimming, while the first one is likely responsible for the dimming itself.

## 5.2. Excitation and damping of higher overtones

As we show in Sect. 3, the higher overtones of the FM period are not observed only after the dimming, but occasionally even before, usually when the photospheric layers are long-term infalling, as in Fig. 5. Furthermore, it appears that full hysteresis loops are not formed when higher overtones dominate. This perhaps suggests that when the material is long-term infalling, different processes from beneath the photosphere are capable of exciting the motion of photospheric material to higher overtones.
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Detail from Fig. 1 showing presence of low-order overtones long before the dimming. The *upper panel* shows variability of brightness and temperature; *lower panel* shows variability of velocity. The second and third red arrows denote the estimated times of maxima, which are delayed only by about 150 days from each other. The first red arrow denotes a possible local maximum outflow velocity of inner layers (C1–C3), while the two outermost layers appear not to have a local maximum of velocity. Thus, the behavior of the atmospheric layers slightly resembles the aftermath of dimming. If we had more data, we could see for how long the low-order overtones remained prominent.</figcaption></figure>

As Joyce et al. (2020) proposed, a non-linear excitation may occur; that is, the layers would eventually accumulate enough energy and momentum to compress the material beneath the photosphere. The compressed material would heat and thus impact convective structure near the surface, altering the energy flow beneath the surface.

Assuming that the strength of the two shocks during the dimming is not unique, the shocks could be the source of additional energy or heat that dissipates into the atmosphere. As Granzer et al. (2022) showed, the mode of pulsations began to change shortly before the dimming, at a similar time to when the first powerful shockwave emerges in our data (Fig. 2), while after the second shockwave the change of pulsations became even more apparent. The shocks should generally appear during a rising part of the light curve (Kravchenko et al. 2019, 2021). The first half of the data (Fig. 1) does not show such a powerful shockwave, although it is not impossible that such a shock may have occurred in a time gap in the data, or before the STELLA data collection began.

**6. Conclusions**

Thanks to the enormous data set of spectral observations by STELLA-SES combined with the tomographic method, we were able to analyze the long-term evolution of different layers of Betelgeuse's photosphere with unprecedented detail. We revealed many insights into the variability of the RSG photospheres, pulsations, formation of shockwaves, its passing through the photosphere, and a correlation of different layers of

photosphere with the brightness variability and properties of the CCFs.

We determined the radial velocity of five different layers of Betelgeuse's photosphere over the past 15 yr using the tomographic method. The phase shift of radial velocities between the neighboring photospheric layers can reach about 10 days (therefore 50 days in total), although during a regular cycle the phase shifts are smaller. The variability of the innermost photospheric layer C1 is approximately in phase with the photometric variations (in filter V); therefore, the phase shifts relative to the brightness are the same. The shock-crossing timescale in the photosphere appears to be longer, up to 75 days, until the shock reaches the outermost layer. The determined relative temperature variability based on a line-depth ratio of V I and Fe I (Gray 2000, 2008) is mostly in phase with photometric variations.

Similarly to Granzer et al. (2022), we find that the photospheric velocity is indeed variable on two primary timescales – the LSP of about 2100 days and FM of about 400 days – while during some cycles the low-order overtones become considerably more significant, especially after the dimming. The low-order overtones seem to be more prominent when the photospheric layers are long-term infalling relative to the CoM of Betelgeuse. This possibly suggests that the mode of pulsations may be excited to higher overtones by processes that take place deeper in the atmosphere.

The timescales of hysteresis loops, which relate the phase shift between variations of velocity and temperature, are similar to the timescale of the FM period, and thus we are able to support that the photometric variations are linked to acoustic waves in a similar manner to that proposed by Kravchenko et al. (2019) for μ Cep. Meanwhile, the oscillations of higher overtones appear to fail to form complete loops, which may be a clue as to why the higher overtones are being damped. Curiously, the LSP also appears to form a hysteresis loop. This could be used as evidence that the LSP period is indeed due to the motion of material in giant convection cells.

We revealed two powerful shockwaves related to the Great Dimming. The first one (reported by Kravchenko et al. 2021) is likely the major progenitor of the dimming, causing a major long-lasting outflow of material (Granzer et al. 2022; Dupree et al. 2022), while the second shockwave emerged after the dimming and contributed to the changes of the atmosphere. The outflow caused by the first shockwave was followed by an extreme infall velocity of all layers. By the time of maximum infall velocity, a second shockwave emerged. This shock was even more powerful than the first one, at least in the outer layers. The second shockwave lead to another significant outflow. These events resulted in major changes to the pulsation characteristics of the photosphere. However, these changes did not take effect simultaneously in all the layers. For about a year and half, the pulsation mode of inner layers was already excited to the first overtone, or to some other low-order overtone, while outer layers remained less affected, apparently continuing their previous movement. The dissonance ended by the beginning of 2022, when all the layers now clearly pulsated with the first overtone. Therefore, the rearrangement of the photosphere appeared complete. Based on the previous behavior of the star, we expect the FM to eventually regain its dominance. Nonetheless, the behavior of the star remains unprecedented, as this was followed by the second strongest expansion of the inner layers (as shown by the gradients) in the beginning of 2023 and a significant brightening of the star in April 2023. Finally, for the last six months (since July, 2023), the star has not exhibited any significant variations of brightness.
The CCF profiles show weak asymmetries rather than clear line doubling, likely related to passing shockwaves or convective motions. Similarly to Kravchenko et al. (2019), we used the FWHMs of CCFs for the detection of shocks in the photoshere. The increased FWHMs revealed two powerful shockwaves during the dimming, which are clearly unprecedented relative to the rest of the data. The FWHMs also show smaller scale variability, which could be related to shocks of lower amplitudes or convective motions.

Our results demonstrate that the triggering mechanism for mass-loss events are powerful shocks, while the higher overtones of radial pulsations may be excited in the process. The episodic gaseous outflows are likely the missing component to explain the mass-loss process in RSGs (Humphreys & Jones 2022).

This study further validates the potential of the tomographic method. Most importantly, it reveals many new major insights into the nature of RSGs. We intend to follow up on these results in a forthcoming work where we will compare the results to 3D radiative-hydrodynamics simulations.

*Acknowledgements.* We thank the editor and anonymous referee for useful comments that improved the quality of the paper. Based on data obtained with the STELLA robotic telescopes in Tenerife, an AIP facility jointly operated by AIP and IAC. We acknowledge with thanks the variable star observations from the AAVSO International Database contributed by observers worldwide and used in this research. D.J. acknowledges support from the ERASMUS+ programme of the European Union and from ESO-MEYS Training Programme of the Czech Republic. K.P. acknowledges funding from the German Leibniz Community under project P67/2018. A.C. acknowledges support from the French National Research Agency (ANR) funded project PEPPER (ANR-20-CE31-0002).

## References

Albrow, M. D., & Cottrell, P. L. 1996, MNRAS, 278, 337

Alvarez, R., Jorissen, A., Plez, B., Gillet, D., & Fokin, A. 2000, A&A, 362, 655

Alvarez, R., Jorissen, A., Plez, B., et al. 2001a, A&A, 379, 288

Alvarez, R., Jorissen, A., Plez, B., et al. 2001b, A&A, 379, 305

Cannon, E., Montargès, M., de Koter, A., et al. 2023, A&A, 675, A46

Chatys, F. W., Bedding, T. R., Murphy, S. J., et al. 2019, MNRAS, 487, 4832

Chiavassa, A., Haubois, X., Young, J. S., et al. 2010, A&A, 515, A12

Chiavassa, A., Freytag, B., Masseron, T., & Plez, B. 2011, A&A, 535, A22

Chiavassa, A., Kravchenko, K., & Goldberg, J. A. 2024, Living Rev. Comput. Astrophys., 10, 2

Davies, B., & Plez, B. 2021, MNRAS, 508, 5757

Dharmawardena, T. E., Mairs, S., Scicluna, P., et al. 2020, ApJ, 897, L9

Drevon, J., Millour, F., Cruzalèbes, P., et al. 2024, MNRAS, 527, L88

Dupree, A. K., Strassmeier, K. G., Matthews, L. D., et al. 2020, ApJ, 899, 68

Dupree, A. K., Strassmeier, K. G., Calderwood, T., et al. 2022, ApJ, 936, 18

Freytag, B., Liljegren, S., & Höfner, S. 2017, A&A, 600, A137

Granzer, T., Weber, M., Strassmeier, K. G., & Dupree, A. 2021, in The 20.5th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun (CS20.5), 41

Granzer, T., Weber, M., Strassmeier, K. G., & Dupree, A. 2022, in Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, 185

Gray, D. F. 2000, ApJ, 532, 487

Gray, D. F. 2008, AJ, 135, 1450

Guinan, E., Wasatonic, R., Calderwood, T., & Carona, D. 2020, The Astronomer's Telegram, 13512, 1

Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951

Harper, G. M., DeWitt, C., Richter, M. J., et al. 2017, ApJ, 836, 22

Harper, G. M., Guinan, E. F., Wasatonic, R., & Ryde, N. 2020, ApJ, 905, 34

Hick, P., Buffington, A., & Jackson, B. V. 2007, SPIE Conf. Ser., 6689, 66890C

Humphreys, R. M., & Jones, T. J. 2022, AJ, 163, 103

Jackson, B., Buffington, A., Hick, P., et al. 2004, Solar Phys., 225, 177

Jadlovský, D., Krtička, J., Paunzen, E., & Štefl, V. 2023, New Astron., 99, 101962

Josselin, E., & Plez, B. 2007, A&A, 469, 671

Joyce, M., Leung, S.-C., Molnár, L., et al. 2020, ApJ, 902, 63

Kervella, P., Decin, L., Richards, A. M. S., et al. 2018, A&A, 609, A67

Kiss, L. L., Szabó, G. M., & Bedding, T. R. 2006, MNRAS, 372, 1721

Kravchenko, K., Van Eck, S., Chiavassa, A., et al. 2018, A&A, 610, A29

Kravchenko, K., Chiavassa, A., Van Eck, S., et al. 2019, A&A, 632, A28

Kravchenko, K., Wittkowski, M., Jorissen, A., et al. 2020, A&A, 642, A235

Kravchenko, K., Jorissen, A., Van Eck, S., et al. 2021, A&A, 650, L17

Liljegren, S., Höfner, S., Eriksson, K., & Nowotny, W. 2017, A&A, 606, A6

López Ariste, A., Georgiev, S., Mathias, P., et al. 2022, A&A, 661, A91

Ma, J.-Z., Chiavassa, A., de Mink, S. E., et al. 2024, ApJ, 962, L36

MacLeod, M., Antoni, A., Huang, C. D., Dupree, A., & Loeb, A. 2023, ApJ, 956, 27

Montargès, M., Cannon, E., Lagadec, E., et al. 2021, Nature, 594, 365

Ogane, Y., Ohshima, O., Taniguchi, D., & Takanashi, N. 2022, Open Eur. J. Variable Stars, 233, 1

Paunzen, E., Supíková, J., Bernhard, K., Hümmerich, S., & Prišegen, M. 2021, MNRAS, 504, 3758

Plez, B. 2012, Astrophysics Source Code Library [record ascl:1205.004]

Ryabchikova, T., Piskunov, N., Kurucz, R. L., et al. 2015, Phys. Scr, 90, 054005

Schwarzschild, M. 1952, in Shock Waves in the Atmospheres of Pulsating Stars, 8 (Cambridge University Press), 811

Stothers, R. B. 2010, ApJ, 725, 1170

Strassmeier, K. G., Granzer, T., Weber, M., et al. 2004, Astron. Nachr., 325, 527

Strassmeier, K. G., Granzer, T., Weber, M., et al. 2010, Adv. Astron., 2010, 970306

Taniguchi, D., Yamazaki, K., & Uno, S. 2022, Nat. Astron., 6, 930

Weber, M., Granzer, T., Strassmeier, K. G., & Woche, M. 2008, SPIE Conf. Ser., 7019, 70190L

Weber, M., Granzer, T., & Strassmeier, K. G. 2012, SPIE Conf. Ser., 8451, 84510K

Wood, P. R. 2000, PASA, 17, 18

Wood, P. R., Olivier, E. A., & Kawaler, S. D. 2004, ApJ, 604, 800
Appendix A: Hysteresis loops

<figure><img src="image_6.png" /></figure>

**Fig. A.1.** Hysteresis loops for photometric cycles. **Upper panel:** Estimated photometric cycles along with the line-depth ratios of V I and Fe I. The ratios are determined from the spectra. Consequently, the ratios also show which part of a photometric cycle is covered by the data. **Lower panels:** Hysteresis loops relate the phase shift between the variations of velocity and temperature. To demonstrate the evolution of material in atmospheric layers more clearly, the values were binned in two-week intervals. Most cycles lack significant parts of the loops due to an insufficient number of measurements (especially in the first half of the data) and observation gaps; therefore, in many cases full loops are not observed. The black arrows indicate the direction of a loop, which almost always is as we would expect, that is, counterclockwise. Continued in Fig. A.2.
<figure><img src="image_7.png" /><figcaption>Fig. A.2. Hysteresis loops for photometric cycles; continued from Fig. A.1.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. A.3. Hysteresis loops for photometric cycles; continued from Fig. A.2.</figcaption></figure>
<figure><img src="image_9.png" /></figure>

**Fig. A.4.** Phase shift between variations of velocity and temperature forms hysteresis loops. The loops are shown for approximate LSP cycles and plotted for different layers of the photosphere. The values were binned in 400 d intervals. The shape of the loop of the first cycle (left panels) is affected by data sampling. Nonetheless, the first cycle also suggests that material is returning to its starting position (except in layer C1).