EVOLUTIONARY TRACKS FOR BETELGEUSE

In spite of its brightness, confusion in the proper motion, variability, asymmetry, and large angular diameter have together made the determination of the distance to this star difficult. A summary of the astrometric data is presented in Table 1. The first parallax distance measurements reported by Hipparcos (∼131 pc) and Tycho (∼54 pc) disagreed by more than a factor of two (ESA 1997). This was well outside the range of quoted errors. However, the more recent VLA-Hipparcos distance to Betelgeuse of ∼197 ± 45 pc (Harper et al. 2008) has been derived from multiwavelength observations. This is the value adopted here as having the greatest accuracy and least distortion from the variability.

Because of the variability of the star during observations, it is difficult to ascribe a mean value and uncertainty to the visual magnitude, luminosity, and temperature. The error bars associated with the quoted apparent visual magnitudes and temperatures in Table 2 are largely a measure of the intrinsic variability of the star. They are therefore not a true measurement error. Hence, to assign an uncertainty to the adopted mean visual magnitudes and temperatures, we simply take the unweighted standard deviation of the various determinations of the mean value. Since the luminosity depends on distance, we simply adopt the luminosity of Harper et al. (2008) based on the revised Hipparcos distance in that paper.

Table 2.  Summary of Luminosity Temperature Data for Betelgeuse

Reference	V	$\mathrm{log}(L/{L}_{\odot })$	Teff (K)
Lee (1970)	0.4	⋯	⋯
Wilson (1976)	0.7	⋯	⋯
Sinnott (1983)	0.5 ± 0.6	${4.67}_{-0.40}^{+0.44}$	⋯
Lambert et al. (1984)	⋯	⋯	3800 ± 100
Cheng et al. (1986)	0.42	⋯	⋯
Gaustad (1986)	⋯	⋯	${3250}_{-120}^{+300}$
Mozurkewich et al. (1991)	0.5	⋯	⋯
Dyck et al. (1992)	0.32 ± 0.24	${4.74}_{-0.26}^{+0.31}$	3520 ± 85
Krisciunas (1992)	0.43 ± 0.14	${4.64}_{-0.27}^{+0.29}$	⋯
Di Benedetto (1993)	⋯	⋯	3620 ± 90
Bester et al. (1996)	⋯	⋯	3075 ± 125
Krisciunas & Luedeke (1996)	0.59 ± 0.24	${4.65}_{-0.23}^{+0.31}$	⋯
Dyck et al. (1996, 1998)	⋯	⋯	3605 ± 43
Wilson et al. (1997)	0.64 ± 0.20	${4.61}_{-0.25}^{+0.29}$	⋯
Perrin et al. (2004)	⋯	4.80 ± 0.19	3641 ± 53
Ryde et al. (2006)	⋯	⋯	3250 ± 200
Harper et al. (2008)	⋯	5.10 ± 0.22	⋯
Ohnaka et al. (2011)	⋯	⋯	3690 ± 54
Adopted	${0.51}_{-0.19}^{+0.13}$	5.10 ± 0.22	3500 ± 200

Download table as:  ASCIITypeset image

Determinations of the angular diameter Θ_disk and associated radius of this star from the observations are summarized in Table 3. This is also complicated by various factors. Red supergiants (RSGs) like α Orionis are extended in radius. As such, their surface gravity is smaller than a main-sequence star like the Sun. This results in extended atmospheres and large convective motion that produces asymmetry along with variability in the surface temperature and luminosity. In addition to the random variability in surface luminosity and temperature and intermittent periodic variability, the diameter of this star seems to change with time (Townes et al. 2009).

Table 3.  Summary of Angular Diameter/Radius Data for Betelgeuse

Ref.	Year Obs.	λ (μm)	${{\rm{\Theta }}}_{\mathrm{disk}}^{{\rm{obs}}}(\mathrm{mas})$	${{\rm{\Theta }}}_{\mathrm{disk}}^{{\rm{corr}}}(\mathrm{mas})$	R/R⊙
Michelson & Pease (1921)	1920	0.575	47.0 ± 4.7	55 ± 7	⋯
Balega et al. (1982)	1978	0.405–0.715	56 ± 11	57 ± 7	${803}_{-234}^{+363}$
	1979	0.575–0.773	56 ± 6	57 ± 7	${803}_{-234}^{+363}$
Cheng et al. (1986)	⋯	⋯	42.1 ± 1.1	⋯	${593}_{-118}^{+182}$
Buscher et al. (1990)	1989	0.633–0.710	57 ± 2	⋯	${802}_{-165}^{+256}$
Mozurkewich et al. (1991)	⋯	⋯	⋯	49.4 ± 0.2	${696}_{-126}^{+195}$
Dyck et al. (1992)	⋯	⋯	⋯	46.1 ± 0.2	${623}_{-113}^{+174}$
Wilson et al. (1992)	1991	⋯	⋯	54 ± 2	${761}_{-158}^{+244}$
Bester et al. (1996)	⋯	⋯	⋯	$56.6\pm 1.0$	⋯
Dyck et al. (1996, 1998)	⋯	⋯	⋯	44.2 ± 0.2	⋯
Burns et al. (1997)	⋯	⋯	⋯	51.1 ± 1.5	⋯
Tuthill et al. (1997)	⋯	⋯	⋯	57 ± 8	${803}_{-267}^{+415}$
Wilson et al. (1997)	⋯	⋯	⋯	58 ± 2	${817}_{-170}^{+260}$
Weiner et al. (2000)	1999	11.150	54.7 ± 0.3	55.2 ± 0.5	⋯
Perrin et al. (2004)	1997	2.200	43.33 ± 0.04	43.76 ± 0.12	620 ± 124
Perrin et al. (2004)	⋯	11.15	55.78 ± 0.04	42.00 ± 0.06	620 ± 124
Haubois et al. (2009)	2005	1.650	44.3 ± 0.1	45.01 ± 0.12	⋯
Neilson et al. (2012)	2005	1.650	⋯	44.93 ± 0.15	955 ± 217
Hernandez-Utrera & Chelli (2009)	2006	2.009–2.198	42.57 ± 0.02	⋯	⋯
Tatebe et al. (2007)	⋯	⋯	48.4 ± 1.4	⋯	⋯
Ohnaka et al. (2009)	2008	2.28–2.31	43.19 ± 0.03	⋯	⋯
Townes et al. (2009)	1993.83	11.150	56.0 ± 1.0	⋯	⋯
	1994.60	11.150	56.0 ± 1.0	⋯	⋯
	1999.875	11.150	54.9 ± 0.3	⋯	⋯
	2000.847	11.150	53.4 ± 0.6	⋯	⋯
	2000.917	11.150	55.8 ± 0.9	⋯	⋯
	2000.973	11.150	54.8 ± 1.0	⋯	⋯
	2001.64	11.150	53.4 ± 0.6	⋯	⋯
	2001.83	11.150	52.9 ± 0.4	⋯	⋯
	2001.97	11.150	52.7 ± 0.7	⋯	⋯
	2006.91	11.150	48.4 ± 1.4	⋯	⋯
	2007.96	11.150	50.0 ± 1.0	⋯	⋯
	2008.09	11.150	49.0 ± 1.5	⋯	⋯
	2008.83	11.150	47.0 ± 2.0	⋯	⋯
	2008.93	11.150	47.0 ± 1.0	⋯	⋯
	2009.05	11.150	48.0 ± 1.0	⋯	⋯
	2009.09	11.150	48.0 ± 2.0	⋯	⋯
Ohnaka et al. (2011)	2009	2.28–2.31	42.05 ± 0.05	42.09 ± 0.06	⋯
Adopted	⋯	⋯	55.64 ± 0.04	41.9 ± 0.06	887 ± 203

Download table as:  ASCIITypeset image

Interferometric measurements can be used to determine an effective uniform disk diameter. However, measurements made at one wavelength must be corrected for the variation in the optical depth with wavelength to yield an effective Rosseland mean radius to compare with the Rosseland radius computed in stellar models.

One must also correct the inferred disk for effects of limb darkening, which tend to diminish the observed radius relative to the Rosseland mean radius. This correction can be of order 10% in visible wavelengths or only as little as ∼1% in the near-infrared (Weiner et al. 2002). The presence of hot spots also tends to diminish the apparent radius (Weiner et al. 2003) by as much as 15%. Moreover, the surrounding circumstellar envelope and dust also complicate the identification of the edge of the star. Furthermore, measurements at different wavelengths lead to different results. Infrared measurements over the past 15 yr even seemed to indicate (Townes et al. 2009) that the radius of this star has been recently systematically decreasing (see, however, Ohnaka 2013).

As a result of all of these complications, one must exercise caution when using measurements of the angular diameter. In Table 3 we quote (when available) the uniform disk Rossland mean radius corrected for limb darkening. Quoted values in the literature are distributed into two distinct groups roughly depending on wavelength. One group λ ∼ 1 μm) is centered around 44 mas (Cheng et al. 1986; Mozurkewich et al. 1991; Dyck et al. 1992; Perrin et al. 2004; Haubois et al. 2006, 2009), while the other ($\lambda \sim 10\quad \mu {\rm{m}}$) is centered around 57 mas (Balega et al. 1982; Buscher et al. 1990; Wilson et al. 1992, 1997; Bester et al. 1996; Burns et al. 1997; Tuthill et al. 1997; Weiner et al. 2000).

Photospheric radii must be corrected for limb darkening and wavelength. These corrections are the smallest for the 11 μm measurements (Weiner et al. 2000). Hence, we adopt a weighted average of the 11.15 μm measurements to obtain 55.6 ± 0.04 mas, which would lead to a photospheric Rosseland mean radius of 56.2 ± 0.04 mas.

However, it is argued quite persuasively in Perrin et al. (2004) that the discrepancy between the lower and higher radii could be accounted for in a unified model that includes the possibility of a warm molecular layer around the star, consistent with that also observed in Mira variables. Applying this correction to the 11.15 μm data leads to a 75% correction (Perrin et al. 2004). This reduces the corrected angular diameter to 41.9 ± 0.04 mas. We adopt this as the best means to deduce a radius to compare with stellar models. This adopted angular diameter, combined with the VLA-Hipparcos distance and uncertainty of 197 ± 45 pc, yields a radius of 887 ± 203 R_⊙. These parameters, along with the observed CNO abundances, lack of s-process abundances (Lundqvist & Wahlgren 2005), and the ejected mass, allow for constraints on stellar models as described below.

The surface elemental and isotopic abundances adopted in this study are summarized in Table 4. Surface elemental C, N, and O abundances for Betelgeuse abundances have been measured by Lambert et al. (1984), while ¹²C/¹³C, ¹⁶O/¹⁷O, and ¹⁶O/¹⁸C isotopic ratios have been reported in Harris & Lambert (1984). As noted in those papers and in the calculations reported here, relative to the Sun, Betelgeuse has an enhanced nitrogen abundance, low carbon abundance, and low ¹²C/¹³C ratio. This is consistent with material that has been mixed to the surface as a result of the first dredge-up phase. As we shall see, this places a constraint on the location of this star on its evolutionary track, i.e., that it must have passed the base of the red giant branch (RGB) and is ascending the RSG phase.

Table 4.  Composition Data for Betelgeuse

Quantity			Ref.
[Fe/H] = +0.1			Lambert et al. (1984)
X⊙	Y⊙	Z⊙
0.71	0.27	0.020	Anders & Grevesse (1989)
X	Y	Z
0.70	0.28	0.024	Corrected for Betelgeuse (see text)
(C)	(N)	(O)
8.41 ± 0.15	8.62 ± 0.15	8.77 ± 0.15	Lambert et al. (1984)
12C/13C	16O/17O	16O/18O
6 ± 1	${525}_{-125}^{+250}$	${700}_{-175}^{+300}$	Harris & Lambert (1984)
Adopted Isotopic Mass Fractions
X(12C)	X(13C)	X(${}^{\mathrm{14,15}}{\rm{N}}$)
1.85 ± 0.80 × 10−3	0.31 ± 0.20 × 10−3	4.1 ± 1.1 × 10−3
X(16O)	X(17O)	X(18O)
6.6 ± 2.0 × 10−3	1.3 ± 0.6 × 10−5	1.1 ± 0.4 × 10−5

Download table as:  ASCIITypeset image

However, to compare the observed abundances with our computed models, some discussion is in order. To begin with, Lambert et al. (1984) report [Fe/H] = +0.1 for this star (where $[X]\equiv \mathrm{log}(X/{X}_{\odot })$). We adopt the Anders & Grevesse (1989) protosolar values ${X}_{\odot },{Y}_{\odot },{Z}_{\odot }=0.71,0.27,0.020$, and assuming that [Fe/H] is representative of metallicity, then $[Z]\;=\;+0.1$. This implies Z = 0.024 for this star. Since the helium mass fraction correlates with metallicity ${\rm{\Delta }}Y/{\rm{\Delta }}Z=3$ (Aver et al. 2013), this implies a helium mass fraction of Y = 0.28, from which one can infer $X=(1-Y-Z)=0.70$. Hence, we adopt this composition. We note, however, that newer photospheric abundances (Asplund et al. 2009) would imply a lower metallicity than that adopted here.

Regarding the C, N, and O abundances, Lambert et al. (1984) reported abundances relative to hydrogen: (C), (N), and (O), where $\epsilon (Z)\equiv \mathrm{log}({\rm{N}}(Z)/{\rm{N}}({\rm{H}})+12$ as given in Table 4. They note that the inferred abundances depend on the measured value of T_eff = 3800 K and the assumed value of $\mathrm{log}g=0.0\pm 0.3$ used in their model atmosphere. Although the temperature is different from the value adopted in the present work, this star shows considerable variability, and this was undoubtedly the appropriate temperature during their observing epoch. Hence, there is no need to correct their abundances for temperature.

Values of ${\epsilon }_{i}$ for individual isotopes are straightforward to evaluate from the isotope ratios given in Table 4. The isotopic mass fractions X_i and estimated uncertainties are then determined from our adopted hydrogen mass fraction X_H, i.e.,

$$\begin{eqnarray}&&\mathrm{log}{X}_{i}={\epsilon }_{i}-12+\mathrm{log}({X}_{{\rm{H}}}A),\end{eqnarray} \tag{ 1 }$$

where A is the atomic mass number.

In Table 5 we summarize variability data concerning mass loss, surface brightness features, and surface turbulent velocities. As in Table 1, the adopted values of these quantities represent a range of possible systematic errors plus the intrinsic variability of the star.

Table 5.  Summary of Mass Loss and Variability Data for Betelgeuse

Ref.	$\dot{M}$ (M⊙ yr−1)	Mej (M⊙ )	(${L}_{\mathrm{spot}}/L-1$)	Tspot −Teff (K)	vc (km s−1)
Knapp & Morris (1985)	6 × 10−7	⋯	⋯	⋯	⋯
Glassgold & Huggins (1986)	4 × 10−6	⋯	⋯	⋯	⋯
Bowers & Knapp (1987)	2.2 × 10−6	⋯	⋯	⋯	⋯
Skinner & Whitmore (1988)	5 × 10−7	⋯	⋯	⋯	⋯
Mauron (1990)	4 × 10−6	⋯	⋯	⋯	3
Marshall et al. (1992)	2 × 10−7	⋯	⋯	⋯	⋯
Young et al. (1993)	5.7 × 10−7	0.042	⋯	⋯	⋯
Huggins et al. (1994)	2. × 10−6	⋯	⋯	⋯	⋯
Mauron & Guilain (1995)	2–4 × 10−6	⋯	⋯	⋯	⋯
Guilain & Mauron (1996)	2 × 10−6	⋯	⋯	⋯	⋯
Noriega-Crespo et al. (1997)	⋯	0.034 ± 0.016	⋯	⋯	⋯
Ryde et al. (1999)	2 × 10−6	⋯	⋯	⋯	⋯
Harper et al. (2001)	3.1 ± 1.3 × 10−6	⋯	⋯	⋯	⋯
Plez & Lambert (2002)	2 × 10−6	⋯	⋯	⋯	⋯
Buscher et al. (1990)	⋯	⋯	10%–15%	⋯	⋯
Wilson et al. (1992)	⋯	⋯	12 ± 3%	⋯	6
Tuthill et al. (1997)	⋯	⋯	11%–23%	⋯	⋯
	⋯	⋯	2%–6%	⋯	⋯
Wilson et al. (1997)	⋯	⋯	13.1%–15.0%	600	⋯
	⋯	⋯	4.0%–7.3%	⋯	⋯
	⋯	⋯	6.1%–8.5%	⋯	⋯
Gilliland & Dupree (1996)	⋯	⋯	⋯	200	⋯
Lobel & Dupree (2001)	⋯	⋯	⋯	⋯	12 ± 1
Lobel (2003a)	⋯	⋯	⋯	⋯	12
Le Bertre et al. (2013)	1.2 × 10−6	0.086	⋯	⋯	14
Adopted	2 ± 1 × 10−6	0.09 ± 0.05	12 ± 12%	400 ± 200	9 ± 6

Download table as:  ASCIITypeset image

Observations and models have deduced (Noriega-Crespo et al. 1997; Ueta et al. 2008; Mohamed et al. 2012) that the combination of mass ejection and the supersonic motion of this star relative to the local interstellar medium has led to the formation of a bow shock pointing along the direction of motion. An estimate of the total ejected mass for this star can be deduced from infrared observations of the surrounding dust.

Noriega-Crespo et al. (1997) first analyzed high-resolution IRAS images at 60 and 100 μm, which indicate the presence of a shell bow shock. They deduced a mass in the shell of

$$\begin{eqnarray}&&{M}_{\mathrm{shell}}=0.042\quad {M}_{\odot }\quad \left(\displaystyle \frac{{F}_{60}}{135\;{\rm{Jy}}}\right)\left(\displaystyle \frac{D}{200\;{\rm{pc}}}\right),\end{eqnarray} \tag{ 2 }$$

where F₆₀ is the measured flux from the shell at 60 μm, which they determine to be 110 ± 21 Jy.

For our adopted distance of 197 ± 45 pc, the mass in this shell is then 0.034 ± 0.016 M_⊙ in the immediate vicinity of Betelgeuse. From this, a current mass-loss rate of around (3 ± 1) × 10⁻⁶M_⊙ yr⁻¹ (Harper et al. 2001) can be deduced.

Moreover, it is now recognized that there are multiple bow shocks (Mackey et al. 2012, 2013) indicating that the mass loss from this star is episodic (Decin et al. 2012) rather than continuous. An additional detached shell of neutral hydrogen has also been detected (Le Bertre et al. 2012) that extends out to 0.24 pc. If one considers this shell, then the total ejected mass increases to 0.086 M_⊙, but the average mass-loss rate reduces to 1.2 × 10⁻⁶M_⊙ yr⁻¹ for the past 8 × 10⁴ yr. For our purpose we will adopt a mass-loss rate of (2 ± 1) × 10⁻⁶M_⊙ yr⁻¹ encompassing both the immediate burst rate of Harper et al. (2001) and the average of Le Bertre et al. (2012).

On the other hand, one cannot be sure how much of the material in the bow shock is interstellar and how much is from the star. Also, one cannot distinguish whether there is more undetected matter from previous mass-ejection episodes, or whether some of the material was ejected while still on the main sequence. Hence, for the total ejected mass, we will adopt a value of 0.09 ± 0.05 M_⊙ with a large conservative estimate of the uncertainty in the total mass loss for this star since it has begun to ascend the RGB. This would correspond to a lifetime of between 0.8 × 10⁵ and 1.4 × 10⁵ yr in the RSG phase. This provides a constraint on the models as discussed below.

We utilize a standard mixing-length theory treatment for convection. The Ledoux criterion (i.e., including chemical inhomogeneities) is utilized to identify convective instabilities, and semiconvection is also included (i.e., mixing in regions that are Schwarzschild unstable though Ledoux stable) as described in the MESA code (Paxton et al. (2013).

In a fully three-dimensional hydrodynamical treatment of convection there can be hydrodynamical mixing instabilities at convective boundaries. This is called convective overshoot. As we shall see below, the outer region of Betelgeuse consists of a rapidly developing convective envelope extending deep into the star. There is also convection in the central helium-burning core. Because of this outer convective envelope, the observed properties of this star can be sensitive to the subtitles of convective overshoot at the boundary of the hydrogen-burning shell. Hence, we have also made a study of the affects of convective overshoot on the models derived here.

In both codes this is accomplished via extra diffusive mixing at the boundaries of convective regions. In the MESA code (Paxton et al. 2011, 2013) convective overshoot is parameterized according to the prescription of Herwig (2000). That is, the MESA code sets an overshoot mixing diffusion coefficient

$$\begin{eqnarray}&&{D}_{\mathrm{OV}}={D}_{\mathrm{MLT}}\mathrm{exp}-2z/({{fH}}_{P}),\end{eqnarray} \tag{ 5 }$$

where D_MLT is a diffusion coefficient derived from mixing-length theory and ${H}_{P}=-P/({dP}/{dr})$ is the pressure scale height. The free parameter f denotes the fraction of the pressure scale height that extends a distance z into the radiative region. A typical value for f for AGB stars is f ≈ 0.015 (Herwig 2000), with a maximum value consistent with observed giants of f < 0.3. We consider this range in the models.

However, as we shall see, although the addition of convective overshoot affects the location of the main-sequence turnoff, it has very little effect on the giant branch. Hence, the observed temperature and luminosity cannot be used to fix the convective overshoot parameter. As we discuss below, however, the observed surface abundances place a constraint on the mixing parameter f. Also as discussed below, the age of this star is slightly decreased when convective overshoot is added. Presumably this is due to changes in thermonuclear burning as material is mixed into radiative zones.

The total ejected mass since arriving at the base of the RGB is a small fraction of the total mass both observationally and in the best-fit models. The mass-loss rate increases considerably as the track approaches the base of the RGB. Therefore, to better identify the location of this star along its track, we can compare the observed accumulated ejected mass along with various integrated mass-loss rates since arriving at the base RGB. We note, however, that the mass loss begins slightly earlier in rotating models (Meynet et al. 2013).

As shown in Table 5, our adopted current observed mass-loss rate is (2 ± 1) × 10⁻⁶M_⊙ yr⁻¹ with a total ejected mass of 0.09 ± 0.05 M_⊙ (Knapp & Morris 1985; Glassgold & Huggins 1986; Bowers & Knapp 1987; Skinner & Whitmore 1987; Mauron 1990; Marshall et al. 1992; Young et al. 1993; Huggins et al. 1994; Mauron & Guilain 1995; Guilain & Mauron 1996; Harper et al. 2001; Plez & Lambert 2002; Ryde et al. 2006; Le Bertre et al. 2012; Humphreys 2013; Richards 2013).

This would correspond to a lifetime of roughly between 3 × 10⁴ and 7 × 10⁴ yr in the RSG phase. However, one expects that the mass-loss rate would have varied during the initial ascent from the base of the RGB. Hence, we have considered the various mass-loss rates given below to integrate the total ejected mass from the stellar models.

Realistically modeling the mass loss from Betelgeuse would be quite complicated. It appears to be episodic (Humphreys 2013) and likely involves coupling with the magnetic field (Thirumalai & Heyl 2012), as well as the normal radiatively driven wind. Indeed, there is evidence of a complex MHD bow shock around this star (Mackey et al. 2012, 2013; Mohamed et al. 2012, 2013). Nevertheless, for comparison of observations with models we have considered various parameterized mass-loss rates (Reimers 1975, 1977; Lamers 1981; Nieuwenhuijzen & de Jager 1990; Feast 1992; Salasnich et al. 1999). For the Reimers (1975) rate,

$$\begin{eqnarray}&&\dot{M}=-4\times {10}^{-13}\eta \displaystyle \frac{L}{{gR}}\quad \quad {M}_{\odot }\;{{\rm{yr}}}^{-1},\end{eqnarray} \tag{ 6 }$$

where L, R, and g are in solar units. The observed rate requires a mass-loss parameter of η = 1.34 ± 1 for a ∼20 M_⊙ star of the adopted L and R. This value is not atypical for giants and is very close to the value inferred by Le Bertre et al. (2012). The Reimers (1977) rate

$$\begin{eqnarray}\mathrm{log}(-\dot{M}) & = & 1.50\mathrm{log}(L/{L}_{\odot })-\mathrm{log}(M/{M}_{\odot })\\ & & -2.00\mathrm{log}({T}_{\mathrm{eff}})-4.74\end{eqnarray} \tag{ 7 }$$

implies a mass-loss rate of 3.32 × 10⁻⁶M_⊙ yr⁻¹, which is again consistent with observed rates and the Reimers (1975) value. On the other hand, the Lamers (1981) rate

$$\begin{eqnarray}\mathrm{log}(-\dot{M}) & = & 1.71\mathrm{log}(L/{L}_{\odot })-0.99\mathrm{log}(M/{M}_{\odot })\\ & & -1.21\mathrm{log}({T}_{\mathrm{eff}})-8.20,\end{eqnarray} \tag{ 8 }$$

gives a present mass loss of 8.80 × 10⁻⁶ M_⊙ yr⁻¹, which is higher by a factor of ∼4 than the adopted current mass ejection rate. The mass-loss rate of de Jager et al. (1988)

$$\begin{eqnarray}\mathrm{log}(-\dot{M}) & = & 1.769\mathrm{log}(L/{L}_{\odot })\\ & & -1.676\mathrm{log}({T}_{\mathrm{eff}})-8.158\end{eqnarray} \tag{ 9 }$$

would predict a mass-loss rate of 8.40 × 10⁻⁶M_⊙ yr⁻¹, which is also higher by a factor of ∼4 relative to our adopted rate. The rate from Salasnich et al. (1999)

$$\begin{eqnarray}&&\mathrm{log}(-\dot{M})=-11.59+1.385\mathrm{log}(L/{L}_{\odot })\end{eqnarray} \tag{ 10 }$$

predicts a mass-loss rate of 2.98 × 10⁻⁵M_⊙ yr⁻¹, which is based on the mass-loss pulsation–period relation of Feast (1992)

$$\begin{eqnarray}&&\mathrm{log}(\dot{M})=1.32\times \mathrm{log}P-8.17.\end{eqnarray} \tag{ 11 }$$

This implies a very high current mass-loss rate of 1.96 × 10⁻⁵M_⊙ yr⁻¹ (using the observed 420-day pulsation period).

Based on this comparison, we adopt the Reimers (1975) rate with $\eta =1.34\pm 1$ as best representing this star. We note, however, that this total mass loss may only be a lower limit to the ejected mass along the RSG phase. This is because not all of the ejected mass may be currently detectable if it was ejected sufficiently far in the past.

Figures 4(a) and (b) show H-R diagrams for the 20 M_⊙ modified EG model and the MESA-code model, respectively. The H-R diagrams are quite similar. The dashed line in panel (b) of Figure 4 shows the H-R diagram for a 20 M_⊙ model calculated with the MESA code with an overshoot parameter f = 0.015. Here one can see that the main effect of the convective overshoot is to increase the luminosity of the main-sequence turnoff and the base of the RGB. It does not, however, affect the luminosity or temperature of the star as it moves up the giant branch. Hence, there is little constraint on convective overshoot from the observed luminosity and temperature. A value of f = 0.3, however, moves luminosity of the base of the RGB all the way up to the observed luminosity. This much overshoot, however, is inconsistent with the observed C and N abundances, as we shall see below.

In both of our best-fit models, the first dredge-up occurs shortly after reaching the base of the RGB. This is when the surface nitrogen is enriched. This means that only models in which the current temperature and luminosity correspond to the ascent up the RSG phase can be consistent with this star.

In spite of the large uncertainties, the adopted mass-loss rate and ejected mass around this star can limit the possible present location of this star along its evolutionary track if we accept that the star has only recently ascended the RSG phase as the best-fit models imply. Hence, we (somewhat arbitrarily) deduce the current age of the star from the amount of mass ejected from the base of the RGB. We checked, however, and the deduced age and properties do not depend on this assumption.

Figure 5 shows a comparison of the total mass ejected as a function of time from the base of the RGB. These curves are based on the Reimers (1975) rate with a mass-loss parameter of η = 1.34 for both the EG (thick solid line) and the MESA (thick dashed line) 20 M_⊙ progenitor models. With our adopted parameters we find that by the time the star reached the main-sequence turnoff it had lost 0.1 M_⊙ in both models, and by the time it reached the base of the RGB it had lost 0.3–0.4 M_⊙. It is possible, however, that our adopted Reimer's mass-loss parameter overestimates the main-sequence mass loss. Winds during the main-sequence phase are very different from the winds during the RSG phase. During the main sequence, winds are mainly radiatively driven, while the RSG winds involve molecules, dust, etc. Hence, these tracks represent upper and lower limits for the mass-loss evolution of this star.

Zoom In Zoom Out Reset image size

The shaded areas indicate the excluded regions based on our adopted uncertainty on the total ejected mass in Table 5. The MESA track indicates a more advanced lifetime as it ascends the RGB. The older age for the MESA models is mostly attributable to the increased opacity at low densities and temperatures in the MESA model. This leads to lower luminosity and hence longer lifetimes. In the context of these models, the amount of mass ejected corresponds to a present age since the zero-age main sequence (ZAMS) of ∼8.0 Myr in the EG model and about 8.5 Myr in the MESA model. For both models the ages are about 0.1 Myr less when convective overshoot is included. Here we adopt the ages in the MESA model without convective overshoot as the most realistic, but we include the EG model results as an illustration of the uncertainty in this age estimate.

3.4.1. Interior

Based on our estimated location in its ascent up the RSG phase, we now examine the interior structure associated with this point in its evolution. Figure 6 summarizes some of the interior thermodynamic properties. The EG and MESA models give almost identical results for the interior thermodynamic properties. In both models the star is characterized at the present time by the presence of a developing carbon–oxygen core up to the bottom of the outer helium core at ∼3–4 M_⊙.

Zoom In Zoom Out Reset image size

The central density has risen to ∼10³ g cm⁻³ with a central temperature of ∼10⁸ K. Outside of the developing C/O core there is a region of steadily decreasing density in the outer helium core and hydrogen-burning shell that extends to the outer envelope consisting of low-density (∼10⁻⁶ g cm⁻³) material. The outer envelope reaches to 90% of the mass coordinate of the star. This is followed by rapidly declining density and the development of an outer surface convective zone.

The best-fit MESA model left the main sequence about 10⁶ yr ago, while for the EG model it was only about 3 × 10⁵ yr ago. Both models reached the base of the RGB about 40,000 yr ago. We followed the star through the final exhaustion of core helium burning in both codes, followed by brief epochs of core carbon, neon, oxygen, and silicon burning until core collapse and supernova at an age of 8.5 Myr since the ZAMS for the MESA code. Our best guess is that the star will supernova in less than ∼100,000 yr (even longer in the EG model). We note, however, that the error ellipse encompasses the entire track, so that the star could be further along in its evolution. The constraint that it has passed the first dredge-up, however, means that the star is ascending the RSG phase. Our result is based on mass loss from the base of the RGB and is therefore a lower limit to how far it has evolved as an RSG.

3.4.2. Composition

Figures 7(a)–(b) show composite plots of isotopic abundances versus interior mass from the best-fit EG and MESA 20 M_⊙ progenitor models. These are compared with the surface isotopic compositions determined by Harris & Lambert (1984) and Lambert et al. (1984). The surface CNO surface abundances computed in both models are quite consistent with the observations as shown by the points in the figures. Hence, the models have correctly evolved the star through the first dredge-up.

Zoom In Zoom Out Reset image size

The present-day interior compositions were similar for the two models except that the interior carbon core is slightly larger in the MESA models owing to the slightly older lifetime. Also, the convective core seems to be less efficiently mixed in the EG model than in the MESA simulation. We attribute this to the more sophisticated mixing treatment in the MESA code (Paxton et al. 2013). Nevertheless, in both models the C/O core has built up to a mass fraction of 40%–50% C and O and extends to about 3 M_⊙. Above this, the He core extends to about 6 M_⊙, while the bottom of the outer convective envelope is at 8 M_⊙ in the EG model but has already descended to ∼10 M_⊙ in the MESA model.

The interior convective properties of the EG and MESA models were quite similar. Figure 8 shows a Kippenhahn diagram of the convective regions over the lifetime of a star for a 20 M_⊙ model with α = 1.8 and a convective overshoot parameter of f = 0.015. An arrow at the bottom indicates our deduced present age for Betelgeuse. Dark shaded regions are unstable to convection by the Ledoux criterion. Lighter regions indicate the range of convective overshoot. One can see from this that the star currently has a rapidly developing outer convective envelope extending deep into the interior. One can also see the recent onset of mass loss for this star.

Zoom In Zoom Out Reset image size

The observed surface abundances, however, can be used to place constraints on convective overshoot at the base of the outer convective envelope. The three panels of Figure 9 illustrate the effects of adding convective overshoot with f = 0.015 to the models. One can see that the main effect of convective overshoot is to extend the bottom of the outer convective envelope from 10 to 12 M_⊙. It also causes the surface abundances of C, N, and O to shift away from agreement, with the biggest discrepancy for ¹⁴N, which changes from agreement to a 3σ discrepancy. The reason for these discrepancies is that convective overshoot tends to minimize the influence of the first dredge-up. Based on this, we conclude that the overshoot parameter for the outer region is constrained to be f < 0.10 at the 2σ confidence level and is most consistent with f = 0.0. Hence, we adopt f = 0.0 for our best-fit models. We note, however, that this constraint is not necessarily valid for the inner convective core, where convection could behave quite differently (e.g., Viallet et al. 2015).

Zoom In Zoom Out Reset image size

It has been suggested (Buscher et al. 1990; Wilson et al. 1997; Freytag et al. 2002; Montargés et al. 2015) that the occurrence of bright spots on the surface of Betelgeuse is the result of convective upwelling material at higher temperature. It is possible to examine whether the occurrence of such features is consistent with simple mixing-length theory in our stellar evolution models. That is, the distance l_c over which a surface convective cell moves is characterized in mixing-length theory by the pressure scale height:

$$\begin{eqnarray}&&{l}_{c}=\alpha \left(\displaystyle \frac{P}{{dP}/{dr}}\right).\end{eqnarray} \tag{ 12 }$$

The condition for which a convective cell reaches the photosphere is

$$\begin{eqnarray}&&{l}_{c}\geqslant R-r,\end{eqnarray} \tag{ 13 }$$

where r is the region from which the convective cell begins its upward motion. The temperature T_hs at which this cell appears on the surface as a hot spot is then given by

$$\begin{eqnarray}&&{T}_{\mathrm{hs}}=T(r)+{\displaystyle \int }_{r}^{R}{\left(\displaystyle \frac{{dT}}{{dr}}\right)}_{\mathrm{ad}}{dr},\end{eqnarray} \tag{ 14 }$$

where the adiabatic temperature gradient is

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{dT}}{{dr}}\right)}_{\mathrm{ad}}=\left(1-\displaystyle \frac{1}{\gamma }\right)\displaystyle \frac{T}{P}\displaystyle \frac{{dP}}{{dr}}.\end{eqnarray} \tag{ 15 }$$

The change in luminosity due to the occurrence of such a spot is then

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}L}{L}={\left(\displaystyle \frac{{l}_{c}}{R}\right)}^{2}{\left(\displaystyle \frac{{T}_{\mathrm{hs}}}{T}\right)}^{4},\end{eqnarray} \tag{ 16 }$$

where we have assumed that l_c also characterizes the size of a convective cell when it reaches the surface. This seems justified by 3D simulations and images (Freytag et al. 2002; Chiavassa et al. 2012; Kervella et al. 2013b) that exhibit convective-cell profiles that are roughly spherical. Indeed, a number of detailed three-dimensional numerical simulations of deep convective envelopes have confirmed the adequacy of mixing-length theory except near the interior boundary layers (Chan & Sofia 1987; Cattaneo et al. 1991; Kim et al. 1996). In particular, they also show approximately constant upward and downward mean velocities (Chan & Sofia 1986), implying a high degree of coherence of the giant cells.

In Buscher et al. (1990) a single bright feature was detected that contributed ∼10%–15% of the total observed flux. In Wilson et al. (1997) the observations were consistent with at least three bright spots contributing a total of 20% to the total luminosity. In their best-fit model the hot spots were taken to have a Gaussian FWHM of 12.5 mas corresponding to as much as one-third of the total surface area or about 10% of the surface per hot spot. Haubois et al. (2009) observed two hot spots attributing to about 10% of the total luminosity. In Freytag et al. (2002), the best fit for their radiation hydrodynamic model has luminosity variations of no more than 30% for the total luminosity, with numerous upwelling hot spots. While their model supports the hot spot theory and most parameters fit the data, they derive a mass of only 5–6 M_⊙, which is inconsistent with the observed luminosity and temperature for this star. Further work by Dorch (2004) using MHD modeling may improve on this inconsistency. Here we point out that a simple mixing-length model is marginally adequate to explain the observed hot spots, as we now describe.

Figure 10 shows the surface convection condition, ${l}_{c}/(R-r)$, as a function of interior radius for a 20 M_⊙ model with convective overshoot and α = 1.8. This shows that the typical size of a convective cell near the surface is less than about 2% of the radius of the star. Hence, the observed surface hot spots that cover as much as ∼10% of the stellar disk would have to result from large fluctuations in the size distribution. The temperature change at the surface from our model is ΔT ≈ 300 K, which for a large fluctuation in spot size would correspond to a 20% change in luminosity and be consistent with the observations.

Zoom In Zoom Out Reset image size

We can also deduce a convective velocity by equating the work done in moving the convective cells to the kinetic energy in the bubbles:

$$\begin{eqnarray}&&{v}_{c}=\left(\displaystyle \frac{\alpha k}{\mu {{\rm{m}}}_{{\rm{H}}}}\right){\left(\displaystyle \frac{T\beta }{g}\right)}^{1/2}{\left[{\rm{\Delta }}\left(\displaystyle \frac{{dT}}{{dr}}\right)\right]}^{3/2},\end{eqnarray} \tag{ 17 }$$

where ${\rm{\Delta }}({dT}/{dr})$ denotes the difference between the temperature gradient in the convective cell and that in the surroundings, α is the mixing-length parameter, and β = 1/2 because we are measuring at the center of the convective cell. We calculate a convective velocity at the surface of ${v}_{c}\sim 12$ km s⁻¹. Lobel & Dupree (2001) deduced a value of ${v}_{c}=9\pm 1$ km s⁻¹. The convective velocity has a dependence on α and can therefore serve as a constraint on α. The best-fit models had α = 1.8–1.9. However, using the value of v_c derived from Lobel & Dupree (2001) and Lobel (2003a, 2003b) would imply that α = 1.4 ± 0.2 near the surface.

An alternate explanation posed by Uitenbroek et al. (1998) suggests that the hot spots result from shock waves caused by pulsations of the stellar envelope. They use Bowen's (1988) density stratification model calculated by Asida & Tuchman (1995). In this model, the stellar envelope pulsates, causing repetitive shock waves that create a density profile that is shallow with respect to the density predicted by hydrostatic equilibrium. Future numerical work should be done to resolve which of the two scenarios is correct for α Orionis.

In addition to the random variability due to upwelling convective hot spots, one also expects Betelgeuse to exhibit regular pulsations. As a massive star ascending the red giant, one expects α Ori to exhibit the regular radial pulsations associated with a long-period variable star. Indeed, Betelgeuse is classified (Samus et al. 2011) as a semiregular variable with an SRC subclassification and a period of 2335 days (6.39 yr). Analyzing the periodicity of this star is difficult, however, as more than one cycle appears to be operating simultaneously and intermittently. In Dupree et al. (1987) and Smith et al. (1989) a pulsational period of 420 days was detected that subsequently disappeared. In Dempsey (2015) a pulsation period of 376 days was deduced using data in the AAVSO database.

It is worthwhile to analyze the implied periodicity from the best-fit model deduced here. Although a detailed model is beyond the scope of the present paper, we can estimate the period to be expected from a one-zone linear adiabatic wave analysis of radial oscillations (Cox & Giuli 1968). For a surface shell on the exterior of a homogeneous star of radius R and mass M and a surface adiabatic equation-of-state index γ, the pulsation period Π is just given by a linearization of surface hydrodynamic equations of motion for a radial perturbation δR,

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}(\delta R)}{{{dt}}^{2}}=-(-3{\gamma }_{\mathrm{ad}}-4)\displaystyle \frac{{GM}}{{R}^{3}}\delta R,\end{eqnarray} \tag{ 18 }$$

which has a solution of simple harmonic motion corresponding to a period of

$$\begin{eqnarray}&&{\rm{\Pi }}=\displaystyle \frac{2\pi }{\sqrt{(3{\gamma }_{\mathrm{ad}}-4)(4/3)\pi G\bar{\rho }}},\end{eqnarray} \tag{ 19 }$$

where $\bar{\rho }$ is the mean density of the star. For our best-fit model the mean density is ${\rho }_{0}={10}^{-7.2}$ g cm⁻³ and the mean adiabatic index is ${\gamma }_{\mathrm{ad}}=1.5$. This implies a pulsation period of ${\rm{\Pi }}\approx 770\;\mathrm{day}$s. This seems reasonably close to the intermittent 420-day pulsation period, particularly given the simplicity of the model and the fact that the observed period may be influenced by higher modes. Clearly, a more detailed nonlinear, nonadiabatic analysis that would also determine the amplitude of the pulsations is warranted.

We attempted a more realistic radial pulsation calculation for our model of α Orionis using both the pulsation code from Hansen & Kawaler (1994) and the one in the MESA code. Although the calculated periods for certain harmonics agree with the observed periodicity of α Orionis, the fundamental mode and first overtone could not be modeled well. This problem may be resolved by considering nonadiabatic pulsations coupled to the convection (Xiong et al. 1998). Results of 3D simulations of RSG stars (Jacobs et al. 1998; Freytag et al. 2002) and Betelgeuse in particular (Chiavassa et al. 2011; Freytag & Chiavassa 2013) confirm the development of large-scale granular convection that generates hot spots. Such convective motion should also drive pulsations and should be analyzed in detail to explain the different pulsation frequencies for α Orionis. Indeed, such large-scale convective motion may also contribute to the secondary periods of this star (Stothers 2010).