A Helium-flash-induced Mixing Event Can Explain the Lithium Abundances of Red Clump Stars

Recently, Kumar et al. (2020) used data from GALAH DR2 (Buder et al. 2018) and Gaia DR2 (Gaia Collaboration et al. 2018) to demonstrate that red clump (RC) stars have an enhancement in ${}^{7}\mathrm{Li}$ that is a factor of $\approx 40$ above the value near the tip of the red giant branch (TRGB). This measurement is enabled by asteroseismic diagnostics that allow the separation of H-shell-burning-only red giant stars and He-core-burning red clump stars (e.g., Bedding et al. 2011; Mosser et al. 2011; Hawkins et al. 2018; Ting et al. 2018). Significant surface lithium enhancement is occurring between the TRGB and the RC, where the notable interior stellar event that occurs in this interval is the helium core flash (Thomas 1967). In this Letter, we demonstrate that mixing in the envelope triggered by the first, strongest He subflash and the action of the Cameron–Fowler process (Cameron 1955; Cameron & Fowler 1971) naturally explains this observation.

We construct stellar evolution models using MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019). We evolve a representative model of a 1 ${M}_{\odot }$ star at solar metallicity.¹ We adopt a mixing length of 1.8 times the pressure scale height. Following Kumar et al. (2020), motivated by past work requiring "extra" mixing in RGB stars (e.g., Charbonnel & Zahn 2007; Angelou et al. 2012), these models include thermohaline mixing using the Kippenhahn et al. (1980) prescription with an efficiency of ${\alpha }_{\mathrm{thm}}=100$. This gives the deep mixing necessary to destroy ${}^{7}\mathrm{Li}$ on the first ascent giant branch. The model is non-rotating, and while rotation may play a role in setting the ${}^{7}\mathrm{Li}$ abundance on the MS and RGB (e.g., Charbonnel & Talon 2005; Charbonnel & Lagarde 2010) and rapid rotation a role in producing Li-rich RGB and RC stars (e.g., Martell et al. 2020), rotation seems an unlikely explanation for a ubiquitous mixing event occurring between the TRGB and RC.

We use the standard MESA pp_and_cno_extras nuclear network, which includes 25 species and the reactions covering the pp-chains and CNO-cycles. The MESA nuclear reaction rates are a combination of rates from NACRE (Angulo et al. 1999) and JINA REACLIB (Cyburt et al. 2010). The key rates ${}^{3}\mathrm{He}{(\alpha ,\gamma )}^{7}\mathrm{Be}$ and ${}^{7}\mathrm{Li}{(p,\alpha )}^{4}\mathrm{He}$ are drawn from NACRE below ${10}^{7}\,{\rm{K}}$ and REACLIB above that temperature. We use the electron-capture rate on ${}^{7}\mathrm{Be}$ from Simonucci et al. (2013), as made available in machine-readable form by Vescovi et al. (2019).² Models are initialized on the pre-main sequence with the Asplund et al. (2009) solar abundance pattern and Z = 0.014. This initializes lithium to the meteoritic abundance $A(\mathrm{Li})=3.26$.

The solid line in Figure 1 shows the standard evolutionary track. This closely matches the MESA result shown as the dashed line in Figure 3 of Kumar et al. (2020). Relative to that track, there is a horizontal shift of a few tenths of a dex in the value of $A(\mathrm{Li})$ as the Kumar et al. (2020) calculation was initialized to match the observed abundance of Li at the main-sequence turnoff. With our adopted mixing during the main sequence, lithium is not depleted to the observed main-sequence level, as some extra mixing at the base of the convection zone is typically invoked to explain the observed values (e.g., Baraffe et al. 2017). This discrepancy is not important for our subsequent discussion of the abundance on the RC as the Li experiences substantial dilution and further destruction during the RGB.

Zoom In Zoom Out Reset image size

It has long been understood (Cameron 1955; Cameron & Fowler 1971) that one way to produce ${}^{7}\mathrm{Li}$ is by having a reservoir of ${}^{7}\mathrm{Be}$, produced via ${}^{3}\mathrm{He}{(\alpha ,\gamma )}^{7}\mathrm{Be}$ at interior temperatures $T\gtrsim {10}^{7}\,{\rm{K}}$, that is mixed out into a cooler convective envelope, where the ${}^{7}\mathrm{Be}$ undergoes an electron capture to form ${}^{7}\mathrm{Li}$. We posit that the He flash initiates some mixing that activates this process.

A detailed graphical summary of a typical helium core flash in a 1 ${M}_{\odot }$ MESA model is shown in Figure 2 of Bildsten et al. (2012). The transition to core He burning occurs through a series of subflashes that occur over a period of $\approx 2\,\mathrm{Myr}$. The first flash is the strongest, with the peak He-burning luminosity exceeding ${10}^{9}{L}_{\odot }$. Subsequent subflashes have more modest peak He-burning luminosities $\approx {10}^{3}\mbox{--}{10}^{4}{L}_{\odot }$. Carrying these high luminosities out from the location of the temperature peak requires convective transport, though because the energy goes into expansion (and hence need not be transported through the entire star), the convective zone remains restricted to the He core.

Internal gravity waves (IGWs) can be excited by turbulent convective motions (e.g., Press 1981; Goldreich & Kumar 1990; Lecoanet & Quataert 2013). Recently, Miller Bertolami et al. (2020) suggested that IGWs excited by the He flash might induce photometric variability in hot subdwarf stars (RG cores without large H envelopes) as they evolve toward the He main sequence (RC equivalent). Here we speculate that these IGWs may induce chemical mixing in the region just above the H-burning shell, possibly either wave-driven mixing or wave-driven heating that triggers convective mixing, linking this region with the convective envelope. We note that the effects of IGWs have been previously invoked to explain lithium evolution in other stellar evolution contexts, including mixing as a cause of Li depletion in F stars (Garcia Lopez & Spruit 1991), and the Li abundances of low-mass, solar-type stars (Montalban 1994; Charbonnel & Talon 2005).

A larger Mach number (${ \mathcal M }$) of the convective motions implies increased excitation of IGWs (Lecoanet & Quataert 2013). During the first subflash, the maximum convective Mach number predicted from mixing length theory can reach ${ \mathcal M }\sim {10}^{-2}$, while the later subflashes only reach Mach numbers of ${ \mathcal M }\sim {10}^{-4}$. Just before the flash, the characteristic luminosity near the base of the H-rich layer, ${L}_{{\rm{H}}}$, is set by the H-burning shell to $\sim {10}^{3.5}\,{L}_{\odot }$ (i.e., the TRGB luminosity). If we take the IGW wave luminosity to be ${L}_{\mathrm{IGW}}\sim { \mathcal M }{L}_{\mathrm{conv}}$, then for this first subflash ${L}_{\mathrm{IGW}}\gg {L}_{{\rm{H}}}$, while subsequent flashes have ${L}_{\mathrm{IGW}}\ll {L}_{{\rm{H}}}$. This motivates our suggestion that the first subflash may be special in its effects.

To demonstrate that such a mixing event can explain the observed Li enhancements, we adopt a simple, ad hoc model of flash-induced mixing. In the H-rich envelope material $(X\gt 0.5)$, we apply a constant diffusive mixing coefficient of ${D}_{\mathrm{mix}}={10}^{12}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$. This diffusivity is additive with any other diffusive mixing processes that are occurring (i.e., convective and thermohaline mixing) and the chosen value is intermediate in magnitude between the two. In order to tie this mixing to the He flash, we activate it only when the helium-burning luminosity is large, specifically when $\mathrm{log}({L}_{\mathrm{He}}/{L}_{\odot })\gt 4$, so that mixing occurs primarily during the first He subflash.

The dotted curve in Figure 1 shows the result of a MESA model including this flash-induced mixing. At the first subflash, the star is enhanced up to $A(\mathrm{Li})\approx 1$ while still near the TRGB luminosity and then experiences only a small amount of depletion from that value as it goes to the RC.

In Figure 2, we show the structure at the location of the H-burning shell just before the core He flash. The total mass of ${}^{7}\mathrm{Be}$ present at this location is $\approx {10}^{-11}{M}_{\odot }$ and is significantly larger than the total mass of ${}^{7}\mathrm{Li}$ in the star. The structure of the H-burning shell is set by the He core properties, which are approximately the same in low-mass ($\lesssim 2{M}_{\odot }$), near-core-flash stellar models, and so the mass of ${}^{7}\mathrm{Be}$ should be approximately the same in all low-mass RGs, independent of stellar mass. If this ${}^{7}\mathrm{Be}$-rich region is mixed, without further synthesis or destruction, into an $\approx 0.5{M}_{\odot }$ convective envelope, then this gives a characteristic abundance $A(\mathrm{Li})\approx 0.6$, approximately the abundance seen at the RC. If the true ${}^{7}\mathrm{Li}$ production process is merely dilution and conversion of this pocket of ${}^{7}\mathrm{Be}$, then since the stellar mass increases with the p-mode large frequency separation, ${\rm{\Delta }}\nu$, the observed surface lithium abundance should decrease with increasing ${\rm{\Delta }}\nu$.

Zoom In Zoom Out Reset image size

Figure 3 shows the evolution of $0.9{M}_{\odot }$ and $1.2{M}_{\odot }$ models. All options except the initial mass are unchanged. Reproducing observations of the lithium abundances on the main sequence and first ascent giant branch would require tuning the mixing during these phases. But despite the differences in lithium abundances before the TRGB, because the structure above the H-burning shell is similar, the flash-induced mixing yields similar abundances on the RC.

Zoom In Zoom Out Reset image size

A characteristic length scale to mix is ${\rm{\Delta }}R\sim {10}^{10}\,\mathrm{cm}$, as the width of the ${}^{7}\mathrm{Be}$-rich region is $\approx 5\times {10}^{9}\,\mathrm{cm}$ and this is located $\approx 5\times {10}^{10}\,\mathrm{cm}$ below the lower boundary of the convective envelope. The mixing must occur on the timescale of the first subflash. In our model ${L}_{\mathrm{He}}\gt {10}^{4}\,{L}_{\odot }$ for $\approx 100\,\mathrm{yr}$, so we take a characteristic duration for the subflash to be ${\rm{\Delta }}t\sim {10}^{9}\,{\rm{s}}$. This suggests a required effective mixing diffusion coefficient ${D}_{\mathrm{mix}}\sim {\left({\rm{\Delta }}R\right)}^{2}/({\rm{\Delta }}t)\sim {10}^{11}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$.

Figure 4 shows tracks for different values of ${D}_{\mathrm{mix}}$ and the threshold luminosity above which this mixing is active. Consistent with the estimate in the previous paragraph, a value of ${D}_{\mathrm{mix}}={10}^{10}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$ leads to a track like the standard case in Figure 1 (i.e., without a flash-induced mixing event). Increasing the strength of the mixing allows for the efficient operation of the Cameron–Fowler process and ${}^{7}\mathrm{Li}$ synthesis well beyond the initial amount of ${}^{7}\mathrm{Be}$, as seen in the track for ${D}_{\mathrm{mix}}\,={10}^{14}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1};{L}_{\mathrm{He}}\gt {10}^{4}\,{L}_{\odot }$. (We note, however, that this particular model then experiences a lithium-depletion event associated with the second He subflash on the way to the RC.) Keeping this same higher mixing value, but decreasing the duration of the mixing by increasing the luminosity threshold to ${L}_{\mathrm{He}}\gt {10}^{6}\,{L}_{\odot }$ (corresponding to a duration of $\approx 4$ yr), results in an intermediate amount of ${}^{7}\mathrm{Li}$ synthesis.

Zoom In Zoom Out Reset image size

The ad hoc diffusivities invoked in the previous section illustrate what is required of any mechanism that operates in rough concert with the first He subflash. Here, we provide some simple, order-of-magnitude estimates to evaluate the potential effectiveness of mixing from IGWs.

Figure 5 shows the internal structure of relevant parts of the star around the time of the peak luminosity of the first He subflash. We discuss each quantity as it enters our estimates, but in all panels the light shaded region indicates where we would impose our ad hoc mixing and the darker shaded region is the location of the ${}^{7}\mathrm{Be}$. As in Figure 2, the mixing region contains a mass of $\approx 2\times {10}^{-3}{M}_{\odot }$.

Zoom In Zoom Out Reset image size

Press (1981) makes the energetic estimate that mixing a region requires overcoming the gravitational potential energy and renewing this energy each thermal diffusion timescale. This corresponds to a specific power $\sim {N}^{2}K$, where N is the Brunt–Väisälä frequency and K is the thermal diffusivity. In the mixing region, this is $\sim {10}^{7}{\rm{erg}}\,{{\rm{s}}}^{-1}\,{{\rm{g}}}^{-1}$, implying that the luminosity required to induce mixing is ${L}_{{\rm{mix}}}\sim {10}^{4}{L}_{\odot }$. This quantity is comparable to the luminosity already transported radiatively in this region (i.e., ${L}_{{\rm{H}}}$). At the time shown in the figure, ${L}_{{\rm{IGW}}}\sim { \mathcal M }{L}_{{\rm{conv}}}\sim {10}^{7}{L}_{\odot }$ and the condition ${L}_{\mathrm{IGW}}\,\gtrsim {L}_{\mathrm{mix}}$ is only satisfied during the first subflash. In the Press (1981) framework, the mass diffusivity associated with wave mixing has an upper limit at the thermal diffusivity. The ad hoc mixing invoked in the previous section is of order the thermal diffusivity in the mixing region, $K\sim {10}^{9}-{10}^{13}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$, implying that the efficient wave mixing is required.

The thermal diffusion time across the mixing region is comparable to the duration of the first subflash. Therefore, an integrated energy argument may be better suited to this case of a time-variable luminosity. The overall IGW energy budget is roughly $\int {dt}\,{ \mathcal M }{L}_{\mathrm{He}}\sim 3\times {10}^{5}\,{L}_{\odot }\,\mathrm{yr}\approx 3\times {10}^{46}\,\mathrm{erg}$, where the integral is completely dominated by the first subflash. The energy required to mix this stably stratified region can be estimated as $\int {{dmH}}_{P}^{2}{N}^{2}\sim {10}^{45}{\rm{erg}}$. Therefore, this mixing is energetically plausible, though it does require tapping a significant fraction of the IGW energy.

Examining the profile of N² in the convective region, we see the characteristic convective frequency is ${\omega }^{2}\sim 3\times {10}^{-5}\,{{\rm{s}}}^{-2}$. We expect these waves to have a characteristic scale of the largest eddies in the convective region (i.e., the pressure scale height), and so these are ${\ell }\approx$few modes with horizontal wavenumber, ${k}_{H}=\sqrt{{\ell }({\ell }+1)}/r\sim 1/r$. The traveling waves can break—and hence lead to mixing as the wave energy is transferred into turbulence—if, as given in Lecoanet et al. (2019),

$$\begin{eqnarray}&&{L}_{{\rm{IGW}}}\gtrsim \displaystyle \frac{4\pi {r}^{2}\rho {\omega }^{4}}{{{Nk}}_{H}^{3}}\sim {10}^{8}{L}_{\odot }.\end{eqnarray} \tag{ 1 }$$

(When we quote a single number estimate of such spatially varying quantities, it is the characteristic value evaluated in the ${}^{7}\mathrm{Be}$ region.) This estimate exceeds the peak value of ${L}_{\mathrm{IGW}}$ implying that the traveling waves do not become nonlinear and break.

When IGW mixing is considered in the Sun, typically ${\omega }^{2}\ll {N}^{2}$ and ${N}^{2}\approx$ constant (e.g., Montalban 1994). The profile of N² shows that neither of these conditions are true in this situation. Waves with the characteristic convective frequency have a turning point in the mixing region. The implies a reduction in the wave luminosity with increasing radius since more of the wave spectrum becomes unable to propagate as gravity waves as N² declines. However, waves that reflect can continue to retain their energy. The thermal damping length for the waves (e.g., Zahn et al. 1997) is

$$\begin{eqnarray}&&{L}_{\mathrm{damp}}=\displaystyle \frac{4\pi {r}^{2}\rho {\omega }^{4}}{{{NN}}_{{\rm{T}}}^{2}K}\sim {10}^{14}\,\mathrm{cm},\end{eqnarray} \tag{ 2 }$$

where N_T is the thermal part of the Brunt–Väisälä frequency. (In the ad hoc mixing region, which is beyond the composition gradient of the H shell, ${N}_{T}\approx N$.) The damping length is much longer than the size of the cavity, so wave energy can accumulate in the cavity and the increased amplitudes associated with standing waves in the cavity may make nonlinear effects more important.

To estimate the diffusivity from mixing by IGWs we follow the approach of Montalban (1994) and Montalban & Schatzman (1996), who find

$$\begin{eqnarray}&&{D}_{M}={A}^{2}{K}^{2}\left(\displaystyle \frac{N}{2\pi }\right){\left(\displaystyle \frac{{\rho }_{0}}{\rho }\right)}^{2}{\left(\displaystyle \frac{{r}_{0}}{r}\right)}^{12}{f}^{-2}{u}_{M}^{5}{\alpha }_{M}^{-1},\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&f={\int }_{{r}_{0}}^{r}{{KN}}^{3}{\left(\displaystyle \frac{{r}_{0}}{r^{\prime} }\right)}^{3}{dr}^{\prime} \end{eqnarray} \tag{ 4 }$$

and quantities with 0 subscripts are evaluated at the convective boundary. We assume the leading numerical constant ${A}^{2}\sim 1$ and adopt a characteristic horizontal eddy scale ${\alpha }_{M}\sim {10}^{9}\,\mathrm{cm}$. We follow their estimate for the characteristic eddy velocity of

$$\begin{eqnarray}&&{u}_{M}\sim {\left(\displaystyle \frac{5{L}_{0}}{4\pi {r}_{0}^{2}{\rho }_{0}}\right)}^{1/3}\sim {10}^{6}\,\mathrm{cm}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 5 }$$

evaluated with ${L}_{0}\sim {10}^{7}{L}_{\odot }$.

We plot the value of Equation (3) in the bottom panel of Figure 5. Because the Montalban (1994) approach assumed ${N}^{2}\gg {\omega }^{2}$, it cannot describe the entire ad hoc mixing region, and we show it only where $N\gt \omega$. At the location of the ${}^{7}\mathrm{Be}$, ${D}_{M}\sim 3\times {10}^{9}\,{{\rm{cm}}}^{2}\,{{\rm{s}}}^{-1}$. This illustrates the plausibility of mixing of the required magnitude from IGWs, but more detailed work specialized to this scenario will be required to evaluate whether IGW-induced mixing can quantitatively explain the observed phenomenon.

We have demonstrated that a brief chemical mixing event in the stellar envelope, triggered at the same time as the first, most vigorous He subflash, can reproduce the observed ${}^{7}\mathrm{Li}$ enhancement of stars on the red clump relative to the first ascent red giant branch. We showed that the amount of ${}^{7}\mathrm{Be}$ present above the H-burning shell just before the flash is sufficient, once mixed into the envelope and converted to ${}^{7}\mathrm{Li}$, to explain the typical abundance on the red clump. Given the large luminosities transported by convection in the core during the helium flash, we suggest that the excitation of IGWs by the turbulent convection can provide a physical mechanism to induce this mixing. Our simple, ad hoc mixing prescription has limited predictive power, but characterizes the required properties of the mixing. The association with the first subflash does predict that objects should be enhanced in lithium for their entire evolution between the TRGB and the red clump. Future work using physical models of wave-induced mixing can determine whether this process can meet the necessary criteria in terms of mixing efficiency and spatial location.

We are grateful to Lars Bildsten, Jim Fuller, Adam Jermyn, and Eliot Quataert for helpful discussions and comments on the manuscript. We thank the referee for suggestions that strengthened the presentation of this material. J.S. is supported by the A.F. Morrison Fellowship in Lick Observatory. We acknowledge use of the lux supercomputer at UC Santa Cruz, funded by NSF MRI grant AST 1828315.

Software: MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019), MESASDK 20.3.1 (Townsend 2020), matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), py_mesa_reader (Wolf & Schwab 2017), MesaScript (Wolf et al. 2017).