The evolutionary stage of Betelgeuse inferred from its pulsation periods

Lobel & Dupree (2000) obtained the surface gravity $\log g=-0.5$ and effective temperature $T_{\rm eff}=3500\pm 100$ K by fitting of synthetic profiles to unblended metal absorption lines in the near-IR spectra. Furthermore, from CO equivalent widths for Betelgeuse Carr et al. (2000) obtained $T_{\rm eff}=3540\pm 260$ K. More recently, Levesque & Massey (2020) obtained $T_{\rm eff}=3600\pm 25$ K from optical spectrophotometry. Based on these spectroscopic $T_{\rm eff}$ determinations for Betelgeuse, we adopt $T_{\rm eff}=3500\pm 200$ K for a constraint for our models, and note that the $T_{\rm eff}$ range is the same as that adopted by Dolan et al. (2016).

As emphasized in Josselin et al. (2000), K-magnitude is more useful in deriving bolometric luminosity of Red supergiants (RSG) rather than V-magnitude, because the K-bolometric correction is insensitive to the effective temperature and surface gravity. In addition, K-magnitude is less affected by interstellar extinction or pulsation. From the 2MASS All-Sky Catalog of Point Sources (Cutri et al., 2003), K-magnitude of Betelgeuse is $-4.378\pm 0.186$ mag. Adopting the distance of $222^{+48}_{-34}$ pc from the new combined radio+Hipparcos astrometric solution obtained by Harper et al. (2017), and K-bolometric correction $2.92\pm 0.16$ (Levesque et al., 2005) for $T_{\rm eff}=3500\pm 200$ K, we obtain $\log L/L_{\odot}=5.18\pm 0.11$ for Betelgeuse (which is similar to $5.10\pm 0.22$ adopted by Dolan et al., 2016).

The semi-regular light variations of Betelgeuse indicate several periodicities ranging from $\sim 2200$ d to $185$ d to be simultaneously excited. Indeed, every period analysis for Betelgeuse light curve yielded more than two periods (Jadlovský et al., 2023; Ogane et al., 2022; Wasatonic, 2022; Granzer et al., 2022; Joyce et al., 2020; Kiss et al., 2006; Stothers & Leung, 1971). Table 1 lists most recent period determinations by various authors. Among five groups of periods, we have selected the four periods $P_{1}$,$P_{2}$,$P_{3}$ and $P_{4}$ as listed in the last line of Table 1. (The 365-day period of Jadlovský et al. (2023) is considered to be one year alias.) Regarding these four periods being caused by radial pulsations of Betelgeuse, we have searched for evolution models which simultaneously excite pulsations of the four periods $P_{1},\ldots P_{4}$ in the range of the error box for Betelgeuse in the HR diagram.

We have applied our non-adiabatic radial pulsation code to selected evolution models which enter the error box on the HR diagram (Figure 1). We have found that the four radial pulsation modes having periods similar to those observed in Betelgeuse (right panel of Fig. 1) are excited in models located in the cooler and luminous part of the error box in the HR diagram as shown in Figure 1. Table 2 lists good models in which all the excited radial pulsation periods agree with the observed four periods. These models are in the carbon-burning (or ending) phase having luminosities in the range of $5.26<\log L/L_{\odot}<5.28$. The mass at ZAMS was $19\,M_{\odot}$ and has been reduced to $\sim 11-12\,M_{\odot}$ at the phase of Betelgeuse.

Note that we need a high luminosity in the upper range of error bar, because the predicted period of the third overtone mode deviates from $P_{4}$ quickly as the luminosity decreases, while $P_{1}$ varies slowly as seen in the right panel of Figure 1. This property of $P_{4}$ sets the lowest luminosity acceptable among models with $\varv_{\rm i}=0.2\varv_{\rm crit}$ at model D (Table 2) which has $\log L/L_{\odot}=5.265$ and $P_{4}=181$ days at the shortest limit of the observed range $185\pm 4$ days. The carbon abundance in the core of model D is 0.171, which will be exhausted in 260 years. At every stage of evolution from model D to B (carbon exhaustion), periods of excited pulsation modes agree with those of Betelgeuse. Model C is chosen as a typical model between B and D.

For the models of $\varv_{\rm i}=0.4\varv_{\rm crit}$ we have listed only Model A in Table 2 although the $P_{4}=178$ days is slightly shorter than the observational limit (i.e., a slightly larger luminosity is needed to perfectly agree with $P_{4}$). The central carbon in model A is almost exhausted, and the luminosity and radius do not change any more to the carbon exhaustion phase.

We note that in the lower luminosity range $\log L/L_{\odot}<2.52$, an additional shorter period mode (fourth overtone) is excited (right panel of Fig. 1), but it is not observed, which additionally supports a high luminosity of Betelgeuse.

The spatial amplitude variation of each excited mode in model C as an example is shown in Figure 2 (top and middle panels) by solid lines as a function of temperature in the envelope. In calculating linear pulsation modes we solve a set of equations with an eigenvalue $\sigma$, expressing the temporal variation of radial displacement as $\delta r\!\exp(i\sigma t)$, in which $\sigma$ and $\delta r$ (sometimes called the eigenfunction) are a complex number and a complex spatial function, respectively. The actual pulsation is expressed as the real part $[\delta r\!\exp(i\sigma t)]_{\rm r}$ so that the imaginary part of $\delta r$ expresses the deviation from the standing wave character in the envelope (i.e., zero points of displacement shift during a cycle of pulsation), and $\sigma_{\rm i}$ (imaginary part of $\sigma$) gives the damping/growth rate of the pulsation mode. In contrast to the nearly adiabatic pulsations of the classical Cepheids, pulsations of supergiants are very non-adiabatic so that the imaginary parts of $\delta r$ and $\sigma$ can be comparable to the real parts.

The top panel of Figure 2 shows that $\delta r$ of the longest period $P_{1}$ ($2181$ d) mode is nearly flat in the outermost layers decreasing gently toward the center, which is exactly the property of radial fundamental mode. For this reason we identify $P_{1}$ as the fundamental mode (rather than LSP). The displacement of $P_{2}$ ($434$ d) mode has a dip at $\log T\approx 4.4$ corresponding to a node in the adiabatic pulsation, which indicates $P_{2}$ to be the first overtone mode. No clear nodes appear in strongly nonadiabatic pulsations, because zero points of real and imaginary parts of $\delta r$ are separated. In the middle panel, $\delta r$ (solid lines) $P_{3}$ and $P_{4}$ change rapidly with a few dips. We identify $P_{3}$ and $P_{4}$ as second and third overtones, respectively. We note that all pulsations are confined in the envelope of $\log T<6$.

Note that $|\delta r|$’s for short period modes (for $P_{3}$ and $P_{4}$) in the middle panel of Figure 2 steeply change near the surface. This indicates that the pulsation energy of these modes leaks at the surface, for which we apply running-wave outer boundary condition in solving the eigenvalue problem. In spite of the leakage of energy, we find that these pulsations still grow (i.e., $\sigma_{\rm i}<0$) because the driving in the inner part exceeds the energy loss at the surface.

A dashed line in the top or middle panel of Figure 2 presents the cumulative work $W$ of each mode. The sign of the gradient $dW/dr$ indicates whether the layer locally drives ($dW/dr>0$) or damps ($dW/dr<0$) the pulsation. Furthermore, whether $W$ at the surface is positive or negative means the pulsation to be excited or damped in the model, respectively. The cumulative work $W(r)$ is calculated as (e.g., Saio & Cox, 1980)

where $P$ and $\rho$ are pressure and density, respectively, while $\delta$ means the pulsational perturbation of the next quantity, and Im$(\ldots)$ and $(^{*})$ mean the imaginary part and the complex conjugate, respectively. Figure 2 indicates that the fundamental mode $P_{1}$ of Betelgeuse is mainly driven in the HeII ionization zone at $\log T\sim 4.5-4.6$, while for other shorter period modes the driving in the H/HeI ionization zone ($4.0\la\log T\la 4.4$) is most important.

The bottom panel of Figure 2 shows the ratio of radiative to total luminosity, $L_{\rm rad}/L_{\rm tot}$, which is considerably smaller than unity in large ranges of the envelope indicating a significant fraction of the total energy flux to be carried by convection. The radiative/convective flux ratio affects the period as well as the driving/damping of the pulsation of a star with high luminosity to mass ratio, $L/M>10^{4}$ (in solar units). Because the thermal time is comparable or shorter than the dynamical time in the envelope of such a luminous star, the pulsational variation of energy flux plays an important role to determine the pulsation period not only to excite or damp the pulsation. In fact, if we obtain pulsation periods of model C (Table,2) neglecting convective flux perturbations (which is sometimes called frozen convection approximation) we obtain periods of 1771, 409, 253, and 192 days. Comparing these periods with periods from the full calculations in model C, 2181, 434, 252, and 184 days, we find including convection flux perturbation to be important only for longer period pulsations in such luminous stars. This indicates that if we fitted the observed $2200$-d period with the fundamental period obtained without convection effects, we would have a $\sim 17\%$ larger radius, because the pulsation period is approximately proportional to $R^{1.5}$. In addition, the fact that the convection effect is serious only in the fundamental mode suggests that even with a pulsation code without convection effects, we would have comparable asteroseismic models by using only overtone modes (i.e., avoiding the fundamental mode).

Figure 3 shows theoretical periodic variations of fractional radius (black), temperature (red), and luminosity (blue) at the surface for the radial fundamental mode $P_{1}$ (2181 d; lower panel) and the first overtone $P_{2}$ (433.6 d; upper panel) excited in model C (Table 2). The amplitudes are normalized as $\delta R/R=1$ at $t=0$ for both cases and fixed by suppressing a possible amplitude growth with time.

The amplitudes of the luminosity and the temperature variations of the $P_{1}$ mode are much smaller than the corresponding amplitudes of the $P_{2}$ mode. In addition, the luminosity maximum of $P_{1}$ mode occurs around the phase of maximum displacement $\delta R/R$, while for $P_{2}$ mode it occurs around minimum $\delta R/R$ when $\delta T/T$ is nearly maximum. The difference comes from the different strength of the non-adiabatic effects, which are stronger in the long-period $P_{1}$ mode than in the $P_{2}$ mode. Strong non-adiabaticity reduces the temperature variation significantly so that the radius effect in the luminosity variation

exceeds the temperature effect. This explains why the luminosity maximum of the $P_{1}$ mode occurs around the maximum displacement $\delta R/R$. In contrast, for the shorter-period $P_{2}$ mode, the temperature perturbation to be large enough to exceed the effect of $\delta R/R$ effect so that the luminosity maximum occurs near the phase of temperature maximum.

All recent period analyses by various authors for the Betelgeuse light curves yielded at least periods $P_{1}$ and $P_{2}$ (see e.g. Table 1). This implies that the main features of the light variations of Betelgeuse should be qualitatively represented by a superposition of the periodic variations of $P_{1}$ and $P_{2}$. In order to confirm the property, we have constructed a synthetic semi-periodic variation by adding the two pulsation modes (Figure 4). We have assumed arbitrarily that the two modes start (at $t=0$) with the same amplitudes to each other. While a linear pulsation analysis assumes the (infinitesimal) amplitude proportional to $\exp(-\sigma_{i}t)$ with the imaginary part of the eigenfrequency $\sigma_{i}$ ($<0$ for an excited mode), we assume a slower growth as $(1-\sigma_{i}t)$ because a rapid growth should be suppressed by non-linear effects for the (observable) finite amplitude pulsation. The radial velocity, RV is calculated as $-(2/3)\delta R/dt$, where a factor of 2/3 presents the projection effect from the spherical surface. The RV variation (top panel) is plotted as red-shift (contraction) upward, and the luminosity variation is plotted in magnitude, $\Delta{\rm mag}\equiv-2.5\log(1+\delta L/L)$, where $\delta R/R$ and $\delta L/L$ are sums of the two pulsation modes.

The $1.2$-yr ($P_{2}$; first overtone) mode variations are modulated by the superposition of the $6.0$-yr ($P_{1}$; fundamental) mode pulsation. The modulation amplitude grows with time as pulsation amplitudes grow. It is particularly interesting that the brightness variation (bottom panel of Fig. 4) has a deep minimum at $t\approx 14.5$ mimicking the Great Dimming of Betelgeuse. To fit the luminosity minimum with the Great Dimming ²²2Although complicated phenomena such as dark spots and dust formation occurred during the Great Dimming as summarized in Wheeler & Chatzopoulos (2023), here for simplicity, we have taken into account only pulsational variations for our qualitative synthetic light curve. we have normalized the luminosity variation such that $\Delta{\rm mag}=1.2$ mag at the minimum ($t=14.5$) by multiplying a factor to $\delta L(t)$. RV(t) and $\delta T_{\rm eff}(t)$ are scaled to be consistent with the luminosity variation. On the synthesized light curve, we have superposed V-band photometry results of Betelgeuse from various sources (indicated in this figure), in which all observed magnitudes are shifted by $-0.5$ mag (assuming the mean magnitude to be $0.5$ mag), and observed time is shifted by $-2005.6$ yr to fit the time of the theoretical minimum (at 14.5 yr) with the Great Dimming at $2020.1$ yr.

During the two cycles of $P_{2}$ prior to the Great Dimming, the synthesized light curve agrees with observation showing minima getting deeper, which seems to indicate the Great Dimming to be caused partially by a constructive interference between the fundamental and the first overtone pulsations (Guinan et al., 2019; Harper et al., 2020). Long before the Great Dimming the low-amplitude modulations associated with the fundamental mode ($P_{1}$) approximately trace the envelope of local maxima of the light variation. However, the phases of observational $P_{2}$ pulsation deviate considerably from the synthetic curve. This might mean that a phase modulation occurred in the first overtone ($P_{2}$) pulsation, or shorter period variations of $P_{3}$ and $P_{4}$, not included in our synthetic light curve, significantly modify the first overtone pulsation.

Just after the Great Dimming, the $P_{2}$ ($434$-d) pulsation suddenly became invisible and is replaced with low-amplitude shorter timescale $\sim 200$-d variations (Dupree et al., 2022; Jadlovský et al., 2023), while the fundamental mode ($P_{1}$) pulsation seems to remain intact. This disappearance of the first-overtone pulsation may suggest that a massive mass loss (Jadlovský et al., 2023; Dupree et al., 2022) at the Great Dimming disturbed significantly the $P_{2}$ pulsation. Since the growth time of the first-overtone pulsation is about 5.4 yr (Fig. 2), the $434$-d pulsation may appear again around 2025. MacLeod et al. (2023) also predicts the re-appearance of the $\sim 400$-d pulsation in a similar timescale (5–10 years) although they regard it as fundamental mode.

The predicted effective temperature variation (middle panel of Figure 4) hardly modulates with the period $P_{1}$. This is understood as that the temperature variation of the fundamental mode (corresponding to $P_{1}$) is suppressed significantly due to the strong nonadiabatic effect as seen in Figure 3. Effective temperature variations of Betelgeuse obtained by Harper et al. (2020) and Taniguchi et al. (2022) are plotted for comparison. While observed results have some wiggles that are probably caused by shorter period pulsations not considered here, the observed range of $T_{\rm eff}$ variations are comparable with our simple two-mode prediction. We note that Wasatonic (2022) and Mittag et al. (2023) also obtained similar variations in the effective temperature of Betelgeuse.

Around the Great Dimming, the synthesized RV curve (top panel) attains maximum (red-shift) 0.1-yr after the minimum brightness. This phase relation agrees with the observed relation (Dupree et al., 2022; MacLeod et al., 2023). Our two-mode synthesis predicts $4.5$ km s^-1 at the maximum, which is comparable with observed values (Kravchenko et al., 2021; Dupree et al., 2022). However, the observed minimum RV (maximum expansion velocity) about $-6$ km s^-1 that occurred $\sim$200-d before the Great Dimming is smaller than our model prediction $-1$ km s^-1. The discrepancy may be consistent with the emergence of a shock in the expansion phase found by Kravchenko et al. (2021) and Dupree et al. (2022), or a breakout of a convective plume as discussed in MacLeod et al. (2023). While various complex phenomena are involved in the RV variations, the diameter of Betelgeuse seems to change by $5\sim 10\%$ on time-scales comparable to pulsation periods (Townes et al., 2009; Taniguchi et al., 2022).

Figure 5 shows the evolutions of luminosity, mass, and CNO abundance ratios at the surface for models with different initial (at ZAMS) rotation velocity given as $\varv_{\rm i}/\varv_{\rm crit}=0.4,0.2,0.1$ ($M_{i}=19\,M_{\odot}$). The second panel from the top shows that the mass-loss occurs mainly in the red-supergiant range (core He burning stage), where up to $~{}6\,M_{\odot}$ is lost. The lost mass is larger for a larger $\varv_{\rm i}$ because of a larger luminosity associated with larger helium core, which was produced by extensive rotational mixing around the convective core during the main-sequence evolution.

Because of the large mass loss in the red-supergiant stage, CNO processed matter emerges to the surface and increases N/C and N/O ratios. The N/C and N/O ratios of the models in the core carbon burning stage, whose pulsation periods agree well with the observed periods of Betelgeuse, are, unfortunately, larger than the observed values obtained by Carr et al. (2000). As a numerical example, the model predicts a N/C that is 0.5 dex larger than measured. This discrepancy seems to indicate that the rotational diffusive mixing is overestimated. The present computation has been performed using specific choices for the diffusion coefficients describing the transport by shear turbulence and meridional currents (see references in Ekström et al., 2012). The choice made in the present computation actually favors the mixing. Other choices would have produced smaller results. As a numerical example, models of 15 M_⊙ with an initial rotation equal to 40% the critical velocity at solar metallicity have been computed with different choices of these diffusion coefficients (see Nandal et al.2023 in preparation). The N/C ratio at the surface when $\log T_{\rm eff}=3.6$ can be reduced up to 0.5 dex changing the expressions of the diffusion coefficients describing the shear mixing. Thus indeed, the solution to that problem can be linked to the specific physics used here for rotational mixing. These models would however present similar structure of their envelope at the red supergiant stage and thus we do not expect this will influence the properties of the pulsation modes during that stage.