Standing on the Shoulders of Giants: New Mass and Distance Estimates for Betelgeuse through Combined Evolutionary, Asteroseismic, and Hydrodynamic Simulations with MESA

Both the ∼400 day and 5.6 yr (2050 day) periods are visible in Figure 1, which shows the longitudinal brightness variations of α Ori over the last 90 yr. These visual brightness estimates were collected in large part by amateur observers and archived by the AAVSO.

Zoom In Zoom Out Reset image size

Examining Figure 1 more closely, we see that the amplitude of the brightness drops corresponding to the ∼400 day pulsation period are about 0.3–0.5 mag in the V band. The difference between these and the 1 mag drop in 2019–20 is clear. We do note, however, that Betelgeuse has undergone other periods of drastic dimming a few times over the last 100 yr. Dimming events of comparable magnitude are visible in Figure 1, for instance, in the mid-to-late 1980s and arguably in the early 1950s. An argument could be made for the existence of a 35–40 yr dimming cycle, particularly if we take into account that the sensitivity of instruments has improved considerably over the last few decades. We note that this 3–4 decade variation is of the same order as the suggested rotational period. While this could potentially be a manifestation of rotational modulation, confirmation will require ongoing observation.

The low amplitude and scarcity of adequate comparison stars make visual estimates less sensitive to smaller changes from one pulsation cycle to another. Digital photometric observations exist for the last three decades, but both the quality and quantity have varied over time. Most of the publicly available data have been archived by the AAVSO and provide good coverage from the mid-1980s to the early 2000s and from 2010 onward. To fill in the gap, we supplement the AAVSO data set with observations taken with SMEI.

2.1.1. SMEI Photometry of Betelgeuse

SMEI was designed to follow coronal mass ejections from the Sun and, in order to do this, stellar signals must be removed from its images. About 6000 stellar sources plus the brightest solar system objects were cataloged and then subtracted from the images. It was soon realized, however, that the source subtraction procedure used by the mission can be processed into time series photometry of the brightest stars in the sky, essentially turning SMEI into one of the early space photometry missions (Buffington et al. 2007; Hick et al. 2007). SMEI observed α Ori from early 2003 to late 2011 with a cadence of 104 min. Each year, data collection was split between the three cameras whose outputs needed to be rectified. Yearly systematics arise from the changing thermal conditions in each of the cameras (Tarrant et al. 2008). Slow degradation of the camera sensitivity is also apparent in the data.

We could not remove the annually repeating instrumental signals directly, as the timescale is on the same order as the variation of α Ori. Therefore, we relied on the ensemble photometry of neighboring stars to derive common instrumental characteristics. We inspected 10 nearby bright stars and selected γ Ori, ε Ori, and 32 Ori to generate a template for the instrumental signals. We calculated a smoothed systematics curve by calculating the medians of the combined relative intensity data of these three stars in 4 day windows placed around every time stamp of α Ori, and for each camera separately.

The rectified SMEI light curve of α Ori is the result of scaling the raw data with the systematics curve and then transforming it to magnitudes using m_SMEI = 10.0 mag as the magnitude zero-point. However, the passband of SMEI is not the V band, therefore requiring that we scale and shift the light curve to match the AAVSO data. To compute the appropriate scaling, we determined the brightness difference for six other stars with M1-2 spectral class in the SMEI catalog to be m_V − m_SMEI ≈ 0.15 mag. We found that we needed to stretch the amplitude by a factor of 1.8 to match the V data points. We then averaged the raw photometry points into 1 day bins. While the shape of the variation could also be passband-dependent to some extent, the scarcity of overlapping V data prohibited us from performing a more detailed comparison. The final light curve is plotted in Figure 2, along with the AAVSO V-band data. The SMEI data confirm that the star did not dim excessively during the LSP minimum occurring in 2007–08.

Zoom In Zoom Out Reset image size

A sample table of the processed and binned SMEI photometry and the scaled V-band values can be found in Table 1. Here, we provide formal errors calculated as the shot noise from the number of electron counts.

Table 1.  Processed SMEI Photometry of α Ori

BJD–	SMEI	SMEI	V	V
2400000 (days)	(mag)	Error	(mag)	Error
52677.959995	0.3759	0.0037	0.5168	0.0037
52678.983194	0.3849	0.0039	0.5330	0.0039
52680.041678	0.3717	0.0043	0.5094	0.0043
52680.959028	0.3801	0.0039	0.5244	0.0039
52681.911649	0.3869	0.0035	0.5365	0.0035
52682.934838	0.3908	0.0036	0.5436	0.0036
52683.887465	0.3991	0.0035	0.5586	0.0035
52684.875377	0.4032	0.0035	0.5662	0.0035
52685.863269	0.4030	0.0035	0.5657	0.0035
52686.851169	0.4068	0.0035	0.5727	0.0035
52687.415625	0.4177	0.0093	0.5929	0.0094
52689.038675	0.4142	0.0042	0.5863	0.0042
...

Note. Observations are corrected for systematics and averaged into 1 day bins. Errors calculated as simple shot noise. V mag is the same light curve, scaled to existing V-band data.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

The photometric light curve reveals a richer set of features than the visual light curve. The SMEI observations, in particular, show both the slow variation from the LSP along with additional smaller, more rapid variations. The SMEI data also put the severity of the recent dimming event in perspective: the brightness of the star did not drop below 1.1 mag in the V band during the last 40 yr, whereas the dip commencing in 2019 November dimmed the star to 1.6 mag in that band. The light curve also highlights some smaller variations on the order of a few hundredths of a magnitude on timescales of days to weeks. Similar variations are present in the SMEI light curves of other nearby stars as well, so we do not consider these to be an intrinsic feature of α Ori.

2.1.2. Frequency Analysis of Observations

We analyzed the frequency spectrum of the photometric light curve with Period04 (Lenz & Breger 2005). We were able to identify the LSP and pulsation frequency regions easily, as shown in Figure 3, but the identification of individual frequency components intrinsic to the star was hindered by the presence of yearly aliases. Most notably, the $-{f}_{\mathrm{LSP}}+1/\mathrm{yr}$ component coincides with the pulsation frequency region. As the FM pulsation mode itself is only slightly longer than one year, its harmonics and/or overtones could coincide with yearly aliases.

Zoom In Zoom Out Reset image size

We first apply a pre-whitening procedure to the data with LSP components. Figure 3 shows that the LSP is not strictly cyclic and that α Ori hovered in a bright state throughout the 2010s. We test combinations of multiple harmonics and subharmonics of the main f_LSP = 0.000423 day⁻¹ frequency (P_LSP = 2365 ± 10 days), which is 15% longer than that determined by Kiss et al. (2006), but still within their uncertainty range. We use the 0.5 and 2.5 f_LSP components for the final fit, which successfully reproduces the deep LSP minima in 1985/1989 and in 2001/2006–7. Non-sinusoidal features in LSP light curves are common for smaller red giants in the Magellanic Clouds: one half of the LSP cycle shows an eclipse-like dip, and the other half resembles a plateau. The model proposed by Soszyński & Udalski (2014) to reproduce this shape assumes a nearby orbiting companion and associated dense cloud. Currently, there is no indication of a companion orbiting α Ori, but some observations suggest that the recent dimming was likely caused by excess dust¹⁰ (Levesque & Massey 2020; Safonov et al. 2020).

We detect two significant frequency components (f_puls1 = 0.002469, f_puls2 = 0.002213 day⁻¹) at the pulsation frequency peak, in agreement with the expected short lifetime of the mode. We likewise detect the first harmonic (2 f_puls) of the stronger pulsation component.

Since the pulsation signal is noncoherent, we fit it with a Lorentzian profile as in Kiss et al. (2006), but in combination with a 1/f component to account for the red noise component of the convective motions (bottom panel of Figure 3). We calculate a pulsation frequency of f_puls = 0.00240 ± 0.00014 day⁻¹ from the peak of the profile, corresponding to a period of P_puls = 416 ± 24 days. We can also use the the full width at half maximum (Γ) of the profile to estimate a mode lifetime of τ = 1/πΓ = 1174 days, or around three pulsation cycles. The mode lifetime matches the value calculated by Kiss et al. (2006); the pulsation period, however, is 7.2% longer, though still within their uncertainty range. We note that Dupree et al. (1987) determined a similar, 420 day, period, but this was based on only three years of photometric observations.

Apparent changes to the period likely arise from (1) the non-coherent nature and short lifetime of the mode and (2) interference with photometric variations caused by convective motions and the evolution of hot spots. Differences of up to 15% among apparent periods calculated from shorter and longer data sets have been found for other RSGs as well (Chatys et al. 2019). Here, we report a new period for the photometric data covering only the last three decades; disentangling the temporal evolution of the pulsation is beyond the scope of this work.

2.1.3. Detection of the First Overtone

A standard means of identifying deviations from an assumed periodic signal is the O − C method.¹² Here we attempt to identify and time the various larger and smaller minima in the light curve. The light-curve data appear to alternate between two states: one defined by deep minima exceeding 0.5 mag (e.g., at JD 49800, 52750, 54000, 54400, 58500, and the dip itself at 58800), the other by shallower and more frequent meandering (e.g., around JD 51500, 53200, 55000 to 57000). However, the annual gaps make it difficult to identify enough minima accurately, and it is thus possible that we simply miss one type during certain intervals. Time spans between consecutive shallow minima can be as short as 60–100 days—much shorter than the FM pulsation period. We see no indication of discrete frequency components corresponding to these intervals in the power spectrum of the star, which suggests they are not high-degree pulsation modes. The timescales and low amplitudes, however, do match the convective turnover times of giant convection cells: our photometric results agree with predictions of timescales from 3D radiative hydrodynamic models and the time-resolved results of spectropolarimetric observations of the surface of the star (Freytag et al. 2002; López Ariste et al. 2018). We therefore attribute these short-timescale, stochastic variations to the photometric effects of convective cell turnover on the surface of α Ori.

The critical observational features of Betelgeuse are summarized in Table 2. We add a final note that observations released since the presentation of the lightcurve in this work show a new brightness minimum occurring less than half a pulsation period later, but at a regular depth of V ∼ 1.0 mag (Dupree et al. 2020; Sigismondi et al. 2020). These observations are consistent with behavior seen in our longitudinal data, particularly those periods in which the short-lifetime mode shifts in phase abruptly. It is likewise consistent with the photometric effects of changing convective cells. Our results thus demonstrate that these two phenomena are sufficient to explain the subsequent, short-term semiregular variability of α Ori without invoking explanations such as multiple opaque dust clouds orbiting the star.

Table 2.  Observational Best Values, Estimated Ranges, and Model-derived Constraints (Where Indicated) for α Ori

Property	Value	Source	Comment
Teff	3600 ± 25 K	Levesque & Massey (2020)	range extended by σtheory to ±200 K
Angular diameter	42.28 ± 0.43	Montargès et al. (2014)	limb-darkened
Radius upper limit	∼1100R⊙	Dolan et al. (2016)	data collated from many sources
Radius lower limit	500R⊙	Dolan et al. (2016)	data collated from many sources
Distance	197 ± 45 pc	Harper et al. (2008)	parallax data adopted by Dolan et al. (2016)
Period of variability	388 ± 30 days	Kiss et al. (2006)	dominant, higher frequency; likely FM
Period of variability	2050 ± 460 days	Kiss et al. (2006)	lower frequency; likely LSP
Period of variability (FM)	416 ± 24 days	this work	SMEI+V data; mode ID from GYRE
Period of variability (O1)	185 ± 14 days	this work	SMEI+V data; mode ID from GYRE
Period ratio O1/FM	0.445 ± 0.041	this work	SMEI+V data
Radius	$764+116,-62{R}_{\odot }$	this work	3σ range; seismic analysis
Initial mass	18–21M⊙	this work	median range; seismic analysis
Present-day mass	16.5–19M⊙	this work	median range; seismic analysis
Distance	168.1 + 27.5, −14.9 pc	this work	3σ range; seismic analysis
Parallax	5.95 + 0.58, −0.85 mas	this work	3σ range; seismic analysis

Note. The temperature constraints reflect the spectroscopically derived temperature from α Ori at its brightness minimum, which is not necessarily reflective of its mean temperature. However, even Levesque & Massey (2020)'s 100 K error bars accounting for decadal variations are more restrictive than the theoretical uncertainty imposed by modeled variations in α_MLT. Though we quote a "best" radius and reference a wide range of values, in practice we do not impose any constraints on the radius when modeling. The range of possible radii derived from the models is smaller than the uncertainties reported by many observers. We quote the initial and present-day mass ranges preferred by our seismic models. Masses considered in the initial grid range from 10 to 26 M_⊙.

Download table as:  ASCIITypeset image

We use GYRE to solve the linearized pulsation equations for high-resolution structural models produced during the RSG phase (Townsend & Teitler 2013). The GYRE program is based on a Magnus Multiple Shooting scheme and provides both adiabatic and nonadiabatic calculations. We consider only adiabatic results in this analysis. Figures 5–8 show results from these calculations.

Zoom In Zoom Out Reset image size

As a track that intersects the observational requirements will typically do so at multiple evolutionary timesteps, we can produce several pulsation profiles per track. Where the models are compatible, we generate synthetic frequency spectra at short intervals. In our frequency modeling, we do not restrict our search to FMs (n_pg = 1,¹⁵ l = 0, m = 0) alone; rather, we note the modes and values of any frequencies in the vicinity of Betelgeuse's two dominant periodicities. However, it is immediately clear from our calculations that periods constituting the primary decaying sequences in both panels of Figure 5 are fundamental mode periods.

To account for the theoretical uncertainty on T_eff, we consider two metrics by which an evolutionary track is compatible with the observations. In the upper panel of Figure 5, we show the periods of adiabatic p-modes versus termination age for a collection of models with a range of initial masses and mixing lengths; here, the structural models used in the seismic analysis have effective temperatures strictly within 3600 ± 25 (Levesque & Massey 2020).

In the lower panel of Figure 5, all models use α_MLT = 2.1, but are checked for consistency against the extended, theoretical temperature constraints described in Figure 4. In both panels of Figure 5, all masses refer to the initial mass of the model and the 416 ± 24 day periodicity is denoted by a blue horizontal band.

In the upper panel of Figure 5, the period–age sequence is tighter and more well-defined, but the results between the upper and lower panels are largely consistent. We note that subsequences comprising stars of particular mass are more apparent in the fixed α_MLT case. This visualization more clearly suggests that, at least for masses in the 15–20 M_⊙ range, there will necessarily be some point along the helium-burning branch during which the star will pass through the appropriate frequency band. However, the temporal window during which this occurs is quite small in the context of evolutionary timescales—on the order of 0.5–1.0 Myr. The requirement that this time frame align with a particular observed temperature ends up being quite restrictive.

Collectively, these results suggest a median, model-derived mass range of 16–21 M_⊙ (with some outliers as high as 24M_⊙), at a resolution of 1 M_⊙. This is broadly consistent with other modelers' results, though our results are more accommodating at the lower-mass end.

We repeat this analysis using as a period constraint the probable overtone observed in our photometry (Section 2.1.2) and confirmed in our synthetic frequency spectra: an O1 mode oscillating with a period of 185 days. In lieu of an independent value, we scale the uncertainty of the fundamental mode to derive an uncertainty on the overtone of 13.5 days, yielding P₁ = 185 ± 13.5. This is shown in the red horizontal bars of Figure 6. Though suggested by our photometric analysis, confirmation of the detection of this mode and its classification required supporting evidence from theoretical spectra. This argument is detailed further in Section 4.2.

Zoom In Zoom Out Reset image size

Figure 7 shows other fundamental parameters as a function of period. The FM and its uncertainties are defined by green vertical bars in all panels. An analog using the O1 period and its uncertainties, defined by dark red vertical lines, is shown in Figure 8.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Models in the upper panels of Figures 5 and 6 span the full range of masses and mixing lengths considered in our grid and additionally vary in the prescribed mass-loss coefficient (η = 0.2–1.0), but they are not pre-restricted by agreement with temperature constraints. Instead, these evolutionary tracks are terminated at arbitrary intervals along the helium-burning branch, with spacing set by the degree of helium exhaustion in the core. This is done to produce a more well-populated sequence that incorporates additional sources of uncertainty in the modeling assumptions. Despite this added theoretical noise, the range of possible radii across all models remains heavily restricted by the observational period constraints.

All models in the lower panels of Figures 5 and 6 intersect the theoretically extended temperature uncertainties (which essentially sets their termination ages) and adopt α_MLT = 2.1 and η = 0.8.

In the uppermost panel of graph 1 in Figure 7 and graph 3 in Figure 8, we show radius as a function of FM and O1 period, respectively. Though there is some scatter in the synthetic data, the radial span of the period-compatible models is very narrow in both cases, especially compared to the range of radial estimates collated in Dolan et al. (2016). We recall from earlier discussion that literature estimates of Betelgeuse's radius range from ∼500 R_⊙ to nearly 1300 R_⊙, whereas the results presented here suggest little possibility outside 700–900 R_⊙. If we interpret the FM period measurements as hard limits, our results suggest a radius for α Ori of 764 R_⊙, with 1σ uncertainties of roughly 30 R_⊙ and non-symmetric limits of R_max = 880 R_⊙ and R_min = 702 R_⊙. Figure 8 suggests consistent but even tighter limits of approximately 700–800 R_⊙ We thus report a 3σ, model-derived radius of ${764}_{-62}^{+116}\,{R}_{\odot }$.

In the middle and lower panels of each graph, we show the models' initial masses and terminal masses, respectively, as a function of period. In these plots, we emphasize those models that also have radii within our 3σ uncertainty bounds with larger, darker markers. Considering all possible observational constraints, we report model-derived estimates for the initial and present-day masses of Betelgeuse as approximately 18–21 M_⊙ and 16.5–19 M_⊙, respectively. Though best mass ranges differ slightly among them, these represent median values collated from the four variants of this test (Figures 7 and 8), encompassing two interpretations of the mixing length–temperature degeneracy for each of P₀, P₁. The lower rows of Table 2 summarize these findings. Taking into account the likelihood of a previous merger event, which would significantly complicate any inferences about the state of Betelgeuse at birth, it is our present-day mass estimates that are most pertinent.

In Section 2.1.2 we presented observational evidence of a new frequency component, f₁, that corresponds to a periodicity of P₁ = 185 days. While strong aliasing caused by annual gaps in the data makes this detection somewhat uncertain, the ratio of this mode to the fundamental suggests that it is the first overtone. Because the literature on seismic models of RSGs is sparse (with most works focusing instead on the low- to medium-mass regime; see Trabucchi et al. 2019), we supplement the photometric analysis with theoretical results from GYRE. To test whether the seismic models agree with the observed period ratio we inspect GYRE's prediction for the frequency of the P₁ (n_p = 2, l = 0, m = 0) mode. With the uncertainties propagated from each period, the relevant range is P₁/P₀ = 0.445 ± 0.041.

We find that GYRE's prediction for P₁ is strongly consistent with 185 days and that the observed period ratio gives mass and age predictions that are self-consistent with the mass and age ranges constrained by the fundamental mode, as demonstrated in Figure 9. When mixing length is varied and temperature is fixed, the age and initial mass ranges inferred from the period ratio are 7–11 Myr and 16–23 M_⊙, respectively.

Zoom In Zoom Out Reset image size

If we consider the progression of the models, we find that the period ratio drops in even younger, higher-mass models, with P₁/P₀ as small as 0.30–0.36 at around 5 Myr. Older, lower-mass models reach the 0.50–0.55 range, in agreement with period ratios observed in semiregular and Mira variables (Kiss et al. 1999; Molnár et al. 2019). We thus conclude that the signal detected in our photometric data is indeed consistent with the first overtone, making Betelgeuse a double-mode star.

Observing multiple modes in a pulsating star can provide stringent constraints on its physical parameters: this is the main principle behind asteroseismology (Aerts 2019). The Fourier analysis in Section 2.1.2 only provided us with a formal error for the single detected frequency peak, but we expect the overtone to have a limited mode lifetime, just like the fundamental mode. As Figures 6 and 8 show, and as is discussed in Section 4.1, the first overtone prefers largely the same models as the fundamental mode.

With the radius of α Ori heavily constrained by the seismic models, we can calculate the distance to the star based on the measured angular diameter. Using an angular diameter of 42.28 ± 0.43 mas and the 3σ uncertainty range of the seismic radius estimate, ${764}_{-62}^{+116}{R}_{\odot }$, we calculate ${168}_{-15}^{+27}$ pc for the distance and $\pi ={6.06}_{-0.85}^{+0.58}$ mas for the parallax. Our values are in agreement with the parallaxes derived entirely or in large part from the Hipparcos measurements (see van Leeuwen 2007 and Harper et al. 2008), and place α Ori nearer to us. It is, however, somewhat in tension with the more recent results based on radio observations, with disagreement at the 1–2σ level (Harper et al. 2017). Figure 10 shows our results in context.

Zoom In Zoom Out Reset image size

This discrepancy could stem from various observational or theoretical shortcomings. One possibility is that the period shift is caused by large-amplitude, nonlinear pulsation. Stellar structure adjusts dynamically to the changes caused by coherent pulsation, which may cause a shift in the eigenfrequencies of the structure. Therefore even if the physical parameters of a linear seismic model agree with those of the star, the calculated and observed periods may not. In the case of p-modes, the radius relates to the pulsation period as R ∼ P^2/3 for a given mass. From this alone, we estimate that the linear period of Betelgeuse should be 500 ± 40 days if its radius is 887 R_⊙, as adopted by Dolan et al. (2016). This means that a ∼20% nonlinear decrease would be required to reproduce the observed 416 ± 24 day pulsation period. Given that we cannot yet evolve RSG hydrodynamic models to full-amplitude pulsation, we cannot infer the period shift directly. For comparison, studies show that pulsating Mira models produce period shifts of up to −23% and +15%, which is of the appropriate order (see, e.g., Lebzelter & Wood 2005; Ireland et al. 2011). However, nonlinear period shifts scale with the pulsation amplitude. As such, the relatively low-amplitude, short-lifetime mode seen in α Ori makes such a large shift implausible. Thus, nonlinear pulsation could be, at best, only partially responsible for the discrepancy we observe.

Another possibility is that the true Rosseland angular diameter of the star is smaller than the 41.8 mas value adopted in the present analysis, and thus the diameter of the photosphere could be as small as 36 mas. However, this would suggest that none of the direct imaging and interferometric observations, including the multi-band models created by Perrin et al. (2004) and Montargès et al. (2014), were capturing the photosphere itself, but, rather, only two separate layers in the extended atmosphere. We consider this very unlikely.

On the other hand, consistency between our parallax estimate and the Hipparcos value could indicate that the astrometric cosmic noise has been underestimated in previous studies, or that it is correlated over the timespan of the astrometric observations.

Harper et al. (2017) illustrates that the addition of a large cosmic noise term to the radio data would be required to bring the combined Hipparcos+radio parallax value close to 6.0 mas. Nevertheless, this could be an indication that various effects, such as photocenter displacements and non-axisymmetric stellar disk shapes caused by convective hot spots, have been previously underestimated. A resolved image obtained from ALMA shows a large hot spot toward the disk limb with a temperature contrast of ΔT ≈ 1000 K, which coincides with the rotational axis but provides no context for the temporal evolution of the photocenter displacements (O'Gorman et al. 2017). Kervella et al. (2018) proposed that this hot spot corresponds to a rogue convection cell that might also be magnetically connected to ongoing mass ejection from the polar region of the star. The presence of one or more persistent spots could mean that the cosmic noise level of the star is higher than reported by Harper et al. (2017), and/or that it needs to be modeled as a correlated noise source. Confirming or ruling out this possibility would require either a sustained, multi-year interferometric observational campaign or the development of accurate 3D physical simulations of supergiant stellar atmospheres that can handle the evolution of hot spots. Finding new ways for Gaia to process heavily saturated stars and thus allowing us to obtain better parallaxes and possibly astrometry of individual hot spots is yet another avenue forward (Sahlmann et al. 2018). However, any of these endeavors would require substantial investment from the respective experts.

Finally, with a seismic parallax in hand, we can estimate a tighter luminosity range based on the Stefan–Boltzmann law. Adopting the strict T_eff uncertainty range reported by Levesque & Massey (2020), we report the luminosity constraints of α Ori to be ${\mathrm{log}}_{10}L={4.94}_{-0.04}^{+0.06}$, where L is in units of L_⊙. The effective temperature range permitted by taking into account theoretical uncertainties extends this range to ${\mathrm{log}}_{10}L={4.94}_{-0.06}^{+0.10}$. We note that if we were to superimpose this range on Figure 4, it would intersect the RSB at the appropriate temperatures for tracks with initial masses between 16 and 21 M_⊙. This range is self-consistent with our other means of estimating mass; however, as mass is a derived quantity, we do not use it to restrict the domain of our seismic grid.

We use MESA (version 8118; Paxton et al. 2015) to solve the stellar structure and evolution equations by means of an implicit hydrodynamics scheme in the Lagrangian formalism. Critically, this means that velocity becomes a dynamical variable similar to density, luminosity, and other physical quantities.

To capture the shock propagation, artificial viscosity q is included in MESA's hydrodynamical scheme (Richtmyer & Morton 1967). However, in our calculations, the global motion of the H-envelope remains subsonic (∼0.5 km s⁻¹) until the end of the simulation. Since this is lower than the typical speed of sound in the atmosphere (∼10 km s⁻¹), capturing shocks is not integral to our work, and hence q acts only as a safeguard to prevent numerical instabilities.

Furthermore, the code uses the energy-conserving, time-discretization scheme of Grott et al. (2005), which ensures that the models at two consecutive timesteps are consistent with each other. We describe the precise configuration for our simulations in more detail in the Appendix, and provide the inlists necessary for reproduction.¹⁶

It is important to note that MESA's implicit hydrodynamics scheme does not include the equations of stellar pulsation directly; rather, MESA has a dedicated module, RSP (Radial Stellar Pulsations; Smolec 2016), for calculating nonlinear pulsating models. However, RSP is not suitable for stars with the luminosity-to-mass ratios (L/M) of giants such as α Ori (Paxton et al. 2019). Previous work has shown that the pulsation modes of stars with L/M ratios similar to that of Betelgeuse have high enough growth rates to induce shocks even if the implicit solver is employed (Heger et al. 1997; Yoon & Cantiello 2010; Paxton et al. 2013; Smolec 2016; Goldberg et al. 2020).

While the main goal of incorporating hydrodynamic simulations into our analysis is to study a canonical model of Betelgeuse's envelope rigorously, non-tailored hydrodynamic simulations also provide a second method of calculating short-order pulsation modes, and thus a means of independently verifying the linear calculations. From a small grid of cursory hydrodynamic runs of varying mass, we can estimate theoretical pulsation cycle lengths directly.

We first conduct an exploratory investigation of the hydrodynamic evolution for a small subset of the models in our grid, restricting it to those with initial masses between 17 and 23 M_⊙. We require that the timestep does not exceed some fixed, small value—typically 5000–10,000 s—and compute the temporal evolution for several decades.

Figure 11 demonstrates the oscillatory behavior of two hydrodynamic models with slightly different initial masses. If not handled correctly, the hydrodynamic models will rapidly expand from their evolutionary initial conditions—in most cases, to nearly double our radial limits—before stable pulsations emerge. This is caused by a discrepancy between the luminosity of the inner boundary of the simulated stellar envelope and the actual stellar luminosity. Thus, over time, the star deposits its energy near the surface, causing it to expand. This can be mitigated by applying relaxation procedures to the initial hydrostatic model (Wood et al. 2004; Nicholls et al. 2009; Ireland et al. 2011; Saio et al. 2015).

Zoom In Zoom Out Reset image size

However, it is still possible to derive the pulsation periods and average radii of these models based upon selected cycles before shocks and/or numerical failure occur, thus providing a cursory but independent validation of the pulsation periods computed with GYRE.

The modeled data are produced by estimating the instantaneous period and radius values from the hydrodynamical models using a combination of quadratic and sine functions fit to short segments of the radial evolution. In this way, we can extract multiple theoretical R, P estimates from one hydrodynamic model. We compare these to a set of direct and inferred observations of pulsation periods and radii of variable stars. The bulk of these data come from the collection of Szatmáry (2004),¹⁷ which contains a variety of variable stars including RR Lyrae, Cepheids, and Miras. In addition, we collate period and radius estimates for a number of other supergiants from the available literature: CE Tau, TV Gem, α Sco, V766 Gem, AH Sco, UY Sct, V602 Car, VY CMa, and KW Sgr (Wasatonic & Guinan 1998; Levesque et al. 2005; Kiss et al. 2006; Ohnaka et al. 2013; Wasatonic et al. 2015; Wittkowski et al. 2017). Of these, pulsation and LSP period measurements are available for α Ori and TV Gem, and KW Sgr has two distinct radii published—these have dashed lines connecting their measurements.

Figure 12 shows a set of synthetic periods and radii derived from the hydrodynamic grid. Also shown in Figure 12 are (1) the P–R sequence constructed from observations of pulsators across a wide mass range (gray dots), (2) the observed FM and LSP periodicities for the small number of red giants listed above (red open circles), and (3) the 416 and 2050-day periodicities of Betelgeuse (red closed circles). Masses of the synthetic stars are indicated via the color bar.

Zoom In Zoom Out Reset image size

It is well known that acoustic modes scale with the average density of the star, which itself largely depends on the radius. We should therefore expect a clear correlation between radius and period, as seen both here and in the linear seismic analysis (Figure 7).

As is clear in the left panel of Figure 12, there is a well-defined P, R sequence spanning RR Lyrae up to synthetic supergiants. The periods and radii extracted from the hydrodynamic models extend the established sequence of pulsating stars to higher radii in a systematic and continuous way, indicating that our models experience p-mode pulsation until our simulations are terminated or the numerics break down. A track of given mass can produce multiple data points by sampling at different times, and hence different degrees of radial expansion, during the hydrodynamic calculation.

In turn, the right panel hints at certain mode classifications for some of the observations. The periods for a number of stars with measured radii fall cleanly on the model sequence, and some fall above: the latter could suggest either pulsations in an overtone or that their radii have been overestimated. Given the complicated circumstellar environment surrounding many supergiants, and our own findings on the radius of α Ori, the latter is a plausible explanation. Finally, the LSP signals are clearly separate from the model sequence, confirming once again that the 2050 day periodicity is not driven by acoustic variations.

Even before more careful modeling of the hydrodynamic evolution, it is evident that pulsation periods emerging naturally in the simulations are of the same order as Betelgeuse's 416 day periodicity. The linear and hydrodynamic seismic analyses both demonstrate that this is α Ori's fundamental mode, a fact which, when combined with other classical observations, places particularly strong constraints on the star's radius.

As stars evolve across the HR diagram, they may undergo mode transitions when a new pulsation mode becomes unstable. At this point, the star can switch to the new mode—a phenomenon observed directly in RR Lyrae stars—or transition to a multimode pulsator.

Mode growth rates have various definitions. In the linear framework, they usually represent the natural timescale of changes in the pulsation energy of the star (Catelan & Smith 2015). Growth rates are sometimes calculated directly from the change in amplitude between successive cycles, but one must keep in mind that in nonlinear calculations, amplitudes do not grow indefinitely. Hence, nonlinear growth rates only agree with the linear values initially, eventually fading back to zero (Yoon & Cantiello 2010). Normalized growth rates are thus scaled with the pulsation frequencies. In the case of, e.g., OGLE–BLG–RRLYR–12245, this mode transition lasted for hundreds of cycles, as is consistent with the small growth rates of the modes (Soszyński et al. 2014). But, in contrast to classical pulsators, semiregular stars have very large growth rates that can lead to strong mode interactions, some of which may even become chaotic (Buchler et al. 1996). As such, it is theoretically possible that Betelgeuse has recently experienced a rapid mode transition, or a rapid increase in amplitude of an overtone mode already present, and that the superposition of the resulting modes created the unusually low brightness minimum seen in 2019 November. It is thus worth investigating whether such a situation can be simulated; however, modeling multimode pulsation in the nonlinear regime is notoriously difficult; at present, this is only reliably reproducible for stars with much lower L/M ratios (Kolláth et al. 2002; Smolec 2016). It is thus beyond the scope of the current paper to investigate such a situation, though this scenario is one we hope to address in a subsequent investigation.

We consider the evolution of a canonical hydrodynamic model whose initial and terminal evolutionary conditions are consistent (as best as possible) with the parameters reported in Section 4. We construct a star with initial mass 21 M_⊙ and terminal mass of 19.54 M_⊙ during core helium burning. In order to force the stabilized radius to be consistent with our reported values, we must inflate the mixing length parameter to α_MLT = 3.0. However, we still require that our hydrodynamic model intersect the theoretical temperature uncertainties described above, 3600 ± 200 K, during its oscillations.

Following the general methodology outlined in Goldberg et al. (2020), we cease tracing the evolution of the innermost 6.0M_⊙, representing the core, at the conclusion of the classical evolutionary run. At this point (often called "core removal" or, more accurately, "core freezing"), a typical timestep in the code is still 100–1000 yr. As a numerical test, we run an evolutionary (nonhydrodynamic) model of a 20M_⊙ star and find that after He exhaustion, it takes approximately 50,000 yr for the luminosity to change by 0.05, whereas our hydrodynamic models show luminosity variations of similar scale over the course of 40 yr. Hence, the core can be considered unchanged to good approximation over the duration of the hydrodynamical simulations. It is thus valid to adopt constant inner boundary conditions set by r_cut such that $M({r}_{\mathrm{cut}})=6.0{M}_{\odot },L=L({r}_{\mathrm{cut}})$ throughout our calculations of the short-timescale, dynamical evolution of the envelope. The value of 6.0M_⊙ is chosen so that 1.0 M_⊙ of the outer He-layer remains along with the entire H-envelope. The He layer, which sets the base of the hydrodynamical simulation, forms a core-envelope structure with the H-envelope, and the higher density core ensures that the oscillation of the envelope does not interact directly with the mass gap.

In order to maintain a stable configuration after ceasing to evolve the core, we allow the model to settle into a hydrostatic approximation before turning on the hydrodynamic solver. To capture the short-timescale motion to adequate resolution, we limit the timestep of the hydrodynamical evolution to a maximum of 10⁴ s. A larger timestep of only ∼10⁵ s can already result in the erasure of modes with a subannual period; this is due to the implicit nature of the hydrodynamical solver. We note that an implicit hydrodynamic scheme is necessary in order to follow the global motion of the star because the relative distances among mass shells near the surface are small. In terms of the Courant timescale, it is ∼10⁶ larger than that required by explicit time discretization (∼0.5 s). Thus, in order to track the motion of the surface with sufficient temporal resolution, implicit hydrodynamics must be employed.

Our canonical model is evolved for a total of ∼10,000 steps from the initiation of hydrodynamics until the point at which stellar expansion begins to disturb the pulsation frequency. We find that beyond ∼30 yr, the motion in the star becomes large enough to interact with the convective layer, causing the timestep to drop as low as 1 × 10² s before the code is unable to evolve the model forward in time. At this point, we stop the simulation.

5.4.1. Global Features

We first discuss the temporal evolution of the critical observables in the canonical model. In the four panels of Figure 13, we plot the luminosity, effective temperature, radius, and surface velocity of the star.

Zoom In Zoom Out Reset image size

The system enters into a state of pulsation with steady but growing cycles a few years after the hydrodynamic solver is switched on. Early in the evolution, quasi-annual oscillatory behavior is present in the luminosity and effective temperature, and the stellar radius exhibits a consistent periodic motion on top of a steady exponential expansion. In this work, we will consider the stellar pulsation only when the motion remains linear; we note that once the behavior becomes nonlinear, the timestep becomes too small to follow the pulsation effectively. Moreover, nonlinear pulsations greatly disturb the profile of the stellar envelope, particularly in terms of opacity and free electron fraction, which makes direct comparison difficult.

In Figure 13, four vertical dotted lines indicate moments at which we compare the instantaneous values of the four quantities. The red and green lines correspond to timesteps where the luminosity is at a local maximum and minimum, respectively. The blue and purple lines correspond to timesteps where the surface velocity is at a local minimum and maximum, respectively.

When the star reaches its brightest point in the pulsation cycle, the effective temperature also reaches its maximum. Concurrently, the star is in its most contracted state (radial minimum) and displays a nearly zero surface velocity. This is consistent with the behavior of a classical harmonic oscillator where the displacement is largest during one cycle. Conversely, the luminosity and effective temperature are minimal when the star is most radially extended, and when the star is maximal in surface velocity, the luminosity, effective temperature, and radius are near their average values. We then observe that as the star continues to expand, a wobble in its motion emerges. As indicated in the hydrodynamic evolutionary tracks of Figure 11, the star will eventually approach a state in which it is capable of developing a hydrodynamical instability in the outer envelope—the exact pulsation properties are themselves a function of evolutionary stage and mass, as shown in Yoon & Cantiello (2010). The evolutionary time at which we launch our hydrodynamic simulations is chosen so as best to reproduce the constraints on the classical parameters found in the previous analysis. We adopt a uniform starting point at which the He mass fraction in the core is 10⁻⁴. Though this is somewhat late in the helium-burning phase and thus does not correspond to the evolutionary phase statistically preferred by our classical models, this starting condition lends itself to more stable hydrodynamic models. We find that when the simulations are initiated at earlier or later burning phases, it is difficult to produce stable, suitable stellar models that fit the radius and luminosity constraints derived in previous sections. As our main priority is to match these parameters, we do not further investigate the effects of evolutionary starting condition. We require only that our hydrodynamical initial conditions yield models that reproduce the radius, luminosity, and mass of the evolutionarily preferred models, as these are the features relevant for studying pressure-driven pulsation in the envelope.

The expansion of the outer radius gradually affects the P–R relation, as the sound speed travel time increases with distance. Analysis of the radius is additionally complicated by (1) how the outermost boundary of the star is defined and (2) radiation pressure outside the photosphere. A rigorous treatment of radiative transport is necessary in the photosphere regime, and so we stop the simulation to analyze the motion only when the radius is below this threshold.

The effective temperature T_eff and luminosity L behave similarly to radius in the simulations. In the former case, this is because MESA calculates the effective temperature directly from the luminosity and radius via ${T}_{\mathrm{eff}}^{4}=(L/4\pi {\sigma }_{{\rm{B}}}{R}^{2})$, where σ_B is the Stefan–Boltzmann constant. Given the slow change of the radius, T_eff primarily mirrors the fast-evolving L.

Regarding the evolution of luminosity, we note that from year 15 onward, the star exhibits periodic motion in its brightness. The early motion is highly regular: as in the preliminary grid of hydrodynamic models (see Figure 11), we observe a quasi-annual rise and fall—the correct timescale for the FM. Near the end of the simulation, the large oscillation begins to trigger non-oscillatory motion in all quantities, and this is responsible for the rapid drop in luminosity.

In Figure 14, we plot the structural profiles of six quantities at points indicated by black circles in Figure 13. The global features of the density and temperature profiles show that the outermost 10% of the stellar mass has a relatively low density, sitting between 10⁻⁹ and 10⁻⁷ g cm⁻³, while the temperature lies between 10^3.5 and 10⁵ K. A small density bump appears at ${\mathrm{log}}_{10}(1-q)={\mathrm{log}}_{10}(1-M(r)/M)=-3$, which is accompanied by a sharp fall in temperature. This occurs in order to maintain hydrostatic equilibrium. We note that a density inversion can occur only when convective mixing is inefficient.

Zoom In Zoom Out Reset image size

It should be noted that a density inversion, in general, cannot exist when there is convection: the mixing will smooth out the density difference. The efficiency of this process depends on the ratio of the timescale of photon diffusion to the dynamical time, as discussed in Jiang et al. (2015). A slow photon diffusion time, as applicable to the stellar interior, implies efficient convection. Hence the density discontinuity is an evanescent phenomenon. However, near the surface where the photon diffusion timescale is much shorter—similar to what we have observed—inefficient convective mixing cannot remove the density inversion. As described in Jiang et al. (2015), this density inversion can form and be destroyed repeatedly in the optically thin regions. Numerical schemes, such as the inclusion of a photon "porosity" by modifying the photon acceleration terms, can suppress this phenomenon in optically thick regions (Paczyński 1969).

When convection is less efficient and has a timescale comparable with the dynamical timescale, the density and pressure gradients can change sign, creating a Rayleigh–Taylor instability. Through mixing, the excess density can gradually reduce via diffusion. However, modeling this phenomenon would require a detailed, time-dependent convective scheme, which is not included in this work. For our case, the density inversion plays a less important role in the luminosity evolution, given that this quantity remains steady from ${\mathrm{log}}_{10}(1-q)=-1$ upward (see the subsequent discussion).

The free electron fraction shows that, up to ${\mathrm{log}}_{10}(1-q)\,=-3$, the matter is partially ionized. Beyond that, the low temperature causes the nuclei to recombine with the free electrons. The opacity profile is richer; rather than falling monotonically like the free electron fraction, we see two major opacity bumps near ${\mathrm{log}}_{10}(1-q)=-1.2$ and −2.5. These correspond to the partial ionization zones of H–He i and He ii, respectively (Cox et al. 1973; Kiriakidis et al. 1992).

The velocity and the luminosity vary dynamically during the pulsation. When the stellar luminosity is mid-phase, the whole envelope is contracting with a constant velocity of ∼0.7 km s⁻¹. The whole star contracts more slowly when it is close to its luminosity maximum or minimum. Meanwhile, the luminosity profiles show that when the star is at its local minimum, the luminosity near the interior of the envelope is lower. The opposite applies during its local maximum. We note further that the motion vanishes approaching the core–envelope interface. This suggests that the observed pulsation is a collective motion of the envelope without interaction with the core.

To further quantify changes in the stellar profile during pulsation, we plot in Figure 15 the ratio of the adiabatic temperature gradients for the profiles at luminosity minimum and medium with respect to luminosity maximum. We observe a sinusoidal-like variation of the profile with about seven nodes, located at ${\mathrm{log}}_{10}q\approx -1,-1.5,-1.9,-2.2,-2.3,\,-2.9$, and −3.5. The star at luminosity maximum shows the highest temperature gradient near the surface at ${\mathrm{log}}_{10}q\approx -2.2$ and −3.1. These correspond to the minima in the temperature ratio profile and the extrema in the opacity ratio profiles.

Zoom In Zoom Out Reset image size

These trends are indicative of the κ-mechanism, and thus explain why the pulsation gradually grows over many periods. In particular, during contraction, the higher opacity prevents the heat from being stored in the deeper layers, which in turn prevents unstable energy extraction by ionization.

From this collection of profiles, we can deduce that the ∼1 yr variation is driven by the collective expansion and contraction of the recombined hydrogen layer. We thus conclude that it is this κ-related interaction driving the fundamental pulsation mode in Betelgeuse.

5.4.2. Literature Comparison

To examine the growth of the oscillation further, we plot the total kinetic energy of the system against time in Figure 16. The kinetic energy, which is dominated mostly by the motion of the atmosphere, is much smaller than both the total energy and the total gravitational energy, which are, on average, −6 × 10⁴⁹ and −1 × 10⁵¹ erg, respectively. These are several orders of magnitude larger than the kinetic energy, which is 10^42–43 erg, in agreement with earlier results (see Cox 1980). We note that it is the outermost q = 0.1–1 that contributes to the atmospheric behavior, including the pulsation. This corresponds to ∼2 × 10⁴⁷ erg when this matter is moving at a speed comparable to the escape velocity.

Zoom In Zoom Out Reset image size

We provide an exponential fit in blue in Figure 16 to characterize the rate of growth of the kinetic energy. A function of ${E}_{\mathrm{kin}}=A\exp ({{bt}}_{\mathrm{year}})$ with parameters A = 1.2 × 10⁴⁰ erg and b = 0.2 provides a good fit to the hydrodynamic component of the simulation. Naively, this suggests that when the oscillation begins to grow, we should expect that the outermost layers of material will be expelled by the pulsations within a time of ∼83 yr. In reality, however, we see this rate level off—see, for example, Cox et al. (1966), who showed this in some early pulsation models. This is because viscous and turbulent dissipation limit the maximum amplitude of the star in the nonlinear regime. Historically, the level of dissipation in 1D pulsation models has been tuned to match the observed amplitudes of RR Lyrae and Cepheid stars and to reproduce double-mode pulsation in the models. This dissipation was first applied via the "artificial viscosity" term and later through the eddy viscosity and other α parameters of time-dependent, turbulent convection (see, e.g., Buchler 1990; Kolláth et al. 1998; Takeuti et al. 1998; Smolec & Moskalik 2008).

We note that the models in Figure 17 stop at about ±2% luminosity variation or less, whereas Betelgeuse itself varies by ±10%–30% in V and about ±20%–30% in near-IR (the latter being more closely representative of Betelgeuse's bolometric variation). As such, there is plenty of time remaining for the kinetic energy and amplitude to grow and eventually saturate at that level, but this is beyond what can be achieved with hydrodynamics today before encountering a numerical runaway episode.

Zoom In Zoom Out Reset image size

We note that when the oscillation becomes strong, heat deposition effects near the surface, where there is a sharp density gradient, become important. The extra heat can change the opacity of the matter by increasing the ionization fraction, resulting in stronger amplification of the pulsation and in turn accelerating the predicted timescale from the first pulsation until mass ejection. Concurrently, however, shock dissipation can convert the kinetic energy of the envelope into heat, while mass loss can also efficiently remove kinetic energy from the star. Heat deposition, shock dissipation, and mass loss interact in complicated ways to form the full dynamical picture of the envelope in the later phases.

Similar ejections of the outer layers in models of luminous Cepheid models, however, are known to be related to numerics rather than physics; see Smolec (2016). Regardless, we do not follow the code until this phase because the timestep becomes prohibitively small (∼100 s). In particular, numerical difficulties arise in the Newton–Raphson iterations, during which the code fails to resolve the formation of convection zones around the shock front. Due to shock compression, the extra heating also invalidates the equilibrium assumptions of the mixing length theory (Vitense 1953). To prevent nonphysical convective behavior from developing near the surface, we set a ceiling for the relative convective velocity with respect to local velocity and the speed of sound in the simulations (see also the Appendix for a related setting in MESA). To limit the steepness of shockwaves and distribute them over multiple zones, explicit pulsation codes like RSP include either artificial viscosity or eddy viscosity terms (or both), but this is only effective up to certain L/M ratios (Stellingwerf 1975; Smolec 2016).

Also at this stage, nonlinear effects become dominant, causing subannual features to appear gradually on top of the linear pulsation. As observations of Betelgeuse do not show periodicities on subannual timescales, we do not consider this phase of pulsation further, though a study of nonlinear pulsation with the dynamical coupling of opacity and ionization will be interesting future work.

By comparing the general features of our hydrodynamical model with the pulsation patterns of Betelgeuse, it becomes clear that the quasi-annual variation is indeed caused solely by the contraction and expansion of the star. It is interesting to note that in this linear oscillation phase, we do not see any evidence of longer-timescale variations, such as the 6 and 35 yr periodicities. In fact, the hydrodynamic simulations never reproduce any of the observed variations besides the currently presented quasi-annual pulsation, even when the initial mass is varied. This implies that these periodicities are driven by some mechanism outside the scope of what 1D hydrodynamic simulations can reproduce, i.e., not the κ-mechanism. It would be interesting to conduct further dynamical studies on how the star relaxes when the opacity effects becomes important; however, work in this domain will require an algorithm to suppress the development of the κ-mechanism so that the pulsation can be sustained without triggering excessive mass loss.

There are similar works in the literature that focus on the pulsational features of massive stellar envelopes using the stellar evolution code described in Langer et al. (1988). In particular, Heger et al. (1997) present the dynamical evolution of massive stars from 10 to 20 M_⊙ and analyze their linear stability. In Figure 18, we plot the phase diagram of our canonical model's ${\mathrm{log}}_{10}(L/{L}_{\odot })$ against ${\mathrm{log}}_{10}{T}_{\mathrm{eff}}$ during the hydrodynamic evolution as a means of comparing directly with Figure 5 in Heger et al. (1997). In their work, the 11 M_⊙ RSG model is followed for about 75 periods of oscillation, whereas ours captures the first 45 periods. Beyond the 45th, our models show numerical instability where the expansion and compression interact with the convective mixing zone, which largely suppresses the timestep and creates nonlinear behavior.

Zoom In Zoom Out Reset image size

Heger et al. (1997) show an approximately circular trajectory that spirals outward from ${\mathrm{log}}_{10}{T}_{\mathrm{eff}}({\rm{K}})\approx 3.52$ and ${\mathrm{log}}_{10}(L/{L}_{\odot })=4.90$, whereas our model shows an elliptical trajectory, vacillating between high L and high T_eff on one side and low L and low T_eff on the other. The outward spiraling in our work and theirs demonstrates that both stars are undergoing dynamical instability with a growing amplitude. As expected, Heger et al. (1997)'s model has a lower period because it is a lower-mass model. This implies a more compact envelope, which allows all 75 periods of oscillations to happen within 30 yr—this is only half the time of our model.

Both models show a clockwise trajectory. Since the radius, temperature, and luminosity are related by the blackbody radiation formula $L=4\pi {\sigma }_{{\rm{B}}}{R}^{2}{T}_{\mathrm{eff}}^{4}$, this means that when the stellar models resume their initial luminosities, the models achieve a higher maximum T_eff (i.e., smaller R) and a lower minimum T_eff (i.e., larger R). These features suggest that T_eff, L, and R achieve their local extrema simultaneously in Heger et al. (1997)'s model, but in our case this relationship is slightly lagged. As shown in Figure 13, our model approaches its local extrema with a non-zero velocity; thus, the stellar radius, which affects T_eff, reaches its local maximum and minimum later than L. Our calculations therefore reproduce the phase lag between the luminosity and the velocity that has been observed in many other, smaller pulsators before (Castor 1968; Szabó et al. 2007).

In Yoon & Cantiello (2010), the hydrodynamical features of a 20 M_⊙ star with a luminosity of ${\mathrm{log}}_{10}(L/{L}_{\odot })=5.05$ and temperature of T_eff = 3198 K are analyzed. Their stellar parameters are similar to ours, where our hydrodynamical model assumes a 21 M_⊙ star with a slightly lower initial hydrostatic luminosity at ${\mathrm{log}}_{10}(L/{L}_{\odot })=5.01$ and T_eff = 4000 K. They model about 50 yr of the stellar pulsation; Figure 2 in their work shows the surface velocity and is comparable to Figure 13 in this work. Approximately 20 pulsation cycles are followed in their work, where a higher period of ∼1000 days is observed. Compared to our ∼400 day period, this indicates that their envelope is more relaxed and expanded. Both Yoon & Cantiello (2010) and our work show a consistent growth of the surface velocity. It takes about five cycles for the surface velocity to reach a 10-fold of amplification, while our model takes much longer—almost 20 cycles. This suggests that the κ-mechanism is less efficient in our model, where the star exhibits behavior closer to adiabatic oscillations than driven oscillations. From the growth of kinetic energy of the system, we can estimate that it takes a further ∼10 yr for the pulsation to grow to a surface velocity comparable with that of Yoon & Cantiello (2010). This would correspond to another 11–15 cycles in our case. Despite this, the robust exponential growth of the pulsation energy (see Figure 16) in both works implies that the pulsation could remove the outermost layers of the H-envelope from the star. However, this is not consistent with observational evidence; the RSGs we have observed pulsate with limited amplitude for several decades. The mass-loss rate of Betelgeuse is observed to be roughly ∼10⁻⁶ M_⊙ yr⁻¹ (Dolan et al. 2016). We may speculate that driven pulsation is partially responsible for this driven mass loss, but the pulsational growth rate is limited by the constant energy dissipation through shock heating and mass ejection. It is also true that the typical pulsation amplitude is not as large as would be expected for a driven wind (Fox & Wood 1982; Wood 2000; Wood et al. 2004).

Whether or not mass loss can be driven depends on the degree of saturation in the pulsation of the surface layers (King et al. 1966). When the mode is permitted to develop, this process can be influential in the formation of circumstellar matter in Type-IIn supernovae for massive stars close to ∼20 M_⊙ (e.g., Smith 2017). However, given the regulated oscillation amplitude observed in a number of RSGs empirically, additional mechanisms not modeled in this work must become dominant in regulating the growth of these oscillation patterns.

The most recent analysis of this kind can be found in Goldberg et al. (2020). The pulsation of an RSG with 16.3 M_⊙ is computed using MESA with the GYRE extension (version 11701). In contrast to the approaches discussed above, their hydrodynamic models involve an initial perturbation to the density distribution to trigger direct pulsation of not only the fundamental mode, but also the first overtone. They obtain a star of ${\mathrm{log}}_{10}(L/{L}_{\odot })=5.2$ and 880R_⊙, which is about 10% larger than our model. As a result, their pulsation shows a fundamental mode with a longer period—about 600 days—and first overtone of ∼300 days. We reiterate that since we do not perturb the density profile at the onset of hydrodynamic evolution, and we assume that all pulsation is triggered by numerical perturbations, it is consistent that we do not see higher-order excited modes alongside the fundamental mode. However, the strength of the quasi-annual variation of Betelgeuse but absence of clear, shorter-timescale periods suggests that we should not be concerned by a lack of overtone activity in our hydrodynamic modeling. The results in Goldberg et al. (2020) demonstrate that an overtone in our model would give rise to significant subannual motion, which is not observed in Betelgeuse. Furthermore, the growth rate depends in part on the L/M ratio, and this factor may not be high enough in our models to excite an overtone.

Goldberg et al. (2020)'s density perturbation approach can also cause the star to pulsate with a much larger amplitude initially, which is not achieved in our work nor in the two analyses presented above (Heger et al. 1997; Yoon & Cantiello 2010). Therefore, it is unclear if their work shows a similar exponential growth as in the literature and in our model.

Thus far, we have presented one model with an initial mass of 21 M_⊙. Now, we consider a series of hydrodynamic models of different initial mass and discuss how the progenitor mass affects the pulsation pattern.

We repeat the simulations by varying the progenitor mass, while fixing the mixing length parameter (see the Appendix for details on the exact configuration) so that we can compare consistently among models. In Table 3, we tabulate the global parameters and pulsation statistics of these models. The data present the following trends: when the progenitor mass increases, the present-day mass M_fin, helium core mass M_He, radius at the end of helium burning R, and its corresponding luminosity L all increase. There is a weak decreasing trend for the effective temperature T_eff. Meanwhile, the time required for the nonlinear pulsation to emerge decreases. We note a severe drop between models of 20.2 and 20.5 M_⊙ during which the associated number of pulsations also decreases sharply. Below M = 19 M_⊙, the oscillation does not amplify significantly within 300 yr, at which point we terminate the simulation.

Table 3.  Global Properties and Pulsation Statistics of the Hydrodynamical Models Studied in this Work

M	α	Mfin	MHe	R	${\mathrm{log}}_{10}L$	Teff	trun	Pulse
18	3	17.12	5.57	550	4.92	4160	>300	N/A
19	3	17.90	6.06	624	5.00	4115	>300	N/A
20	3	18.80	6.47	655	5.04	4117	166.2	∼230
20.2	3	18.95	6.60	659	5.05	4120	141.8	∼200
20.5	3	19.17	6.75	707	5.10	4081	31.5	43
21	3	19.54	7.00	721	5.11	4083	27.0	40
22	3	20.30	7.46	787	5.18	4053	21.0	24
23	3	20.92	8.04	875	5.25	4008	16.7	18
19	2.5	17.78	6.07	724	5.01	3832	102.8	122
20	2.5	18.57	6.59	794	5.08	3801	41.1	37

Note. M, M_fin, and M_He are the initial, final, and He-core masses of the star in units of M_⊙, respectively. R, ${\mathrm{log}}_{10}L$, and T_eff are the initial radius in units of R_⊙, luminosity in units of L_⊙ and the effective temperature in units of K at the beginning of the hydrodynamical phase, respectively. α is the mixing length parameter. t_run is the time the pulsation of the star becomes nonlinear, where we stop the simulations, in years. "Pulse" is the number of pulsation cycles experienced by the star before the onset of nonlinear pulsations. No number is available when t_run is larger than 300 yr. All hydrodynamic models are launched from the evolutionary point at which the helium mass fraction in the core is 10⁻⁴. Though this does not correspond to the evolutionary phase statistically preferred by our classical models, this starting condition lends itself to more stable hydrodynamic models and allows us to explore the appropriate radius, luminosity, and mass regimes.

Download table as:  ASCIITypeset image

We note in particular two entries in Table 3 showing simulations with initial (evolutionary) masses of 20M_⊙: one with α_MLT = 3.0 and one with α_MLT = 2.5. In the latter case, the nonlinear excitation time drops considerably. This data point disrupts an otherwise monotonically decreasing trend in pulse number with increasing present-day mass (above 19M_⊙), demonstrating that the relationship between pulse number, t_run, and present-day mass is not totally straightforward. It is well known, however, that below a certain level of precision, the impact on global parameters caused by changes in the mixing length are indistinguishable from variations in mass and metallicity in models of low- to intermediate-mass stars with convective envelopes (Joyce & Chaboyer 2018a). Further, the late-stage evolution of high-mass stars is especially sensitive to convective parameters; in practice, α_MLT is often tuned arbitrarily until the model converges or behaves as desired. We include the last row of Table 3 to highlight this degeneracy and caution against over-interpretation.

In Figure 17, we present the time evolution of the pulsation pattern for models with the progenitor mass from 19 to 22 M_⊙. We choose these masses as their timescales are more relevant to that of Betelgeuse. Before nonlinearity disturbs the simulation, all models behave similarly in both luminosity and radius.

Despite the fact that the excitation time apparently depends on the progenitor mass and mixing length parameter, the the means by which the star becomes excited—i.e., the pulsation driving mechanism itself—is less sensitive to these choices. The dynamical pulsation always concludes with a significant drop in the stellar luminosity, and the peak luminosity and maximum radius are similar among all models near the end of the simulation.

To further outline the similarity, we plot in Figure 18 the phase diagram of representative models from 19 to 23 M_⊙. Clear similarity can be seen for models above 19 M_⊙. In particular, for M = 20 M_⊙, the model has a highly extended t_run of ∼166 yr. All models have an elliptical structure, which is actually a clockwise outward-going spiral. They show once again that all stars evolve toward a high-L and a high-T_eff state simultaneously, or the converse. This suggests that the driving mechanism in all of these models is qualitatively the same, too. A higher progenitor mass gives rise to a sparser trajectory; however, we notice that for M = 19 M_⊙, there is no regularity in the trajectory. This suggests that the κ-mechanism fails to stimulate residual numerical noise into periodic motions.

To further characterize the runaway time of the M = 20 M_⊙ model, we compare the amount of time needed for the star to develop nonlinear pulsation (t_run) after using the hydrodynamical prescription. For progenitor masses above 20.5 M_⊙, t_run decreases slowly with time. As shown in Figure 17, nonlinear activation timescales for masses of this range are between 15 and 40 yr. However, below 20.5 M_⊙, there is a sudden jump in t_run, and the star requires more than roughly ≥150 in order for nonlinearity to become significant. The sudden jump could signify some qualitative changes in the stellar profile, namely that the κ-mechanism becomes much less effective in amplifying the acoustic wave inside the star; a detailed comparison to and analysis of the means of formation for the κ-mechanism will be an interesting future project, but is beyond the scope of the present work. Crucially, this mass-sensitive timescale bifurcation suggests that the time required for the star to develop nonlinear pulsation could be a highly discerning attribute among models of Betelgeuse.

As an order of magnitude estimation, the typical luminosity of our model star is 10⁵ L_⊙. The amount of energy dissipated is then ∼10⁴⁶ erg per year, but the kinetic energy is only on the order of 10^41–43 erg. This is because radiation acts as a damping force through photon emission, and without a consistent driving force for the pulsation, the oscillation would quickly dissipate.

We demonstrated above that there are multiple periodicities in Betelgeuse's lightcurve. These include a quasi-annual mode, a 6 yr mode, a 30 yr modulation and, potentially, an overtone mode with 185 day period. In the hydrodynamic models, we recover only the 416 day period. These results are largely self-consistent, as the 416 day mode is driven by the κ-mechanism, the LSP is not, and the 30 yr modulation is most likely caused by rotation, which is not an internally driven form of variability. In the case of the overtone, however, we must address the question of how multiple modes may appear in the first place.

One possibility is by nonlinear mode excitation, as touched upon in Section 5.3. Through large-amplitude oscillations, the outer layers can accumulate sufficient energy and momentum to compress matter beneath the stellar surface. This results in compression heating, which in turn raises the local temperature. This may impact the convective structure in the near-surface regions, thus presenting an additional source of energy that alters the net energy flow inside the star. Capturing this scenario numerically is particularly challenging because it involves modeling the dynamics of mixing behavior in the convection zone. Meanwhile, the standard mixing length theory adopted in our work assumes the convective mixing is in equilibrium (Vitense 1953). Modeling this phase properly would require a more sophisticated approach to time-dependent mixing and a robust solving mechanism. We reiterate that the development of the overtone mode can also be sensitive to the numerical setting. In particular, the implicit nature of the hydrodynamics tends to suppress acoustic waves naturally, regardless of the use of artificial viscosity. Therefore, the current non-detection can be attributed to numerics alone. Future exploration using, for example, an explicit hydrodynamical scheme would shed light on this matter but is beyond the scope of this study.

We observe that in our hydrodynamical grid, models with M > 21 M_⊙ develop nonlinear pulsation much more quickly than the lower-mass models, as shown in Figure 17. We emphasize that our results only suggest that the star is developing nonlinear pulsation and/or near-surface shocks; how strong the final shock is remains unclear. Because we do not follow our hydrodynamical simulations far enough to investigate rigorously the various energy dissipation channels discussed in Section 5.4, it remains to be studied whether α Ori's strong pulsations eventually reach a new equilibrium or behave in some other way.

Another possible excitation mechanism is wave-driven pulsation, as described in Shiode & Quataert (2014), Fuller (2017), and Fuller & Ro (2018). This mechanism proposes that a convective wave in the star can partially penetrate the evanescent regions¹⁸ and approach the stellar surface. Although wave-driven pulsation was described in the context of very late phases of stellar evolution in those works (i.e., neon–oxygen burning, rather than helium), the theory suggests that as long as convection is activated, energy can be transferred from the interior convection zone to regions near the surface, where it can then excite surface motion. However, depending on the convective luminosity, such a mechanism would provide a heavily condensed energy deposition near the surface, in turn triggering enormous losses in mass of 0.01–1 M_⊙ yr⁻¹. Mass loss of this order is not observed in Betelgeuse.