Is Thermohaline Mixing the Full Story? Evidence for Separate Mixing Events near the Red Giant Branch Bump

The abundances of mixing--sensitive elements including lithium, [C/N], and 12C/13C are known to change near the red giant branch bump. The explanation most often offered for these alterations is double diffusive thermohaline mixing in the stellar interior. In this analysis, we investigate the ability of thermohaline mixing to explain the observed timing of these chemical depletion events. Recent observational measurements of lithium and [C/N] show that the abundance of lithium decreases before the abundance of [C/N], whereas numerical simulations of the propagation of the thermohaline mixing region computed with MESA show that the synthetic abundances drop simultaneously. We therefore conclude that thermohaline mixing alone cannot explain the distinct events of lithium depletion and [C/N] depletion, as the simultaneity predicted by simulations is not consistent with the observation of separate drops. We thus invite more sophisticated theoretical explanations for the observed temporal separation of these chemical depletion episodes as well as more extensive observational explorations across a range of masses and metallicities.


INTRODUCTION
As low-mass stars evolve into red giants, the surface abundances of some elements change in ways that depend on the internal stellar structure. On the early giant branch, the surface convection zone is deepening, diluting the abundance of delicate elements like lithium and beryllium that only survived the main sequence in a narrow region near the surface of the star (Sweigart et al. 1989;Charbonnel et al. 1994;Boothroyd & Sackmann 1999;Boesgaard et al. 2020). As the surface convection zone deepens to its greatest extent in mass, it penetrates into regions that have previously undergone nuclear processing and cycles the processed material to-wards the surface during the first dredge up, leaving behind a chemical discontinuity (see, e.g., Iben 1964;Lambert 1981;Joyce & Chaboyer 2015 and references therein). Due to the nuclear fusion rates in the carbon, nitrogen, oxygen (CNO) cycle and the temperature dependence of this depth, this process will alter the [C/N] and 12 C/ 13 C ratios in a mass-dependent way (Iben 1967;Gratton et al. 2000;Masseron & Gilmore 2015;Martig et al. 2016;Ness et al. 2016). In standard stellar models without additional mixing or physics, the abundances at the surface of the star do not change between this point and the end of the red giant branch (Chanamé et al. 2005;Palacios et al. 2006).
After the first dredge-up, the surface convection zone recedes in mass coordinate as the hydrogen burning shell advances. Eventually the hydrogen burning shell crosses the residual chemical discontinuity and encounters the mean molecular weight gradient left behind by the deepest extent of the surface convection zone (Kippenhahn & Weigert 1994;Fusi Pecci et al. 1990;Charbonnel et al. 1994;Joyce & Chaboyer 2015). At this point, the struc-ture readjusts, the core contracts and the luminosity decreases. In observed stellar populations, this is seen as an over-density of stars and identified as the red giant branch bump (RGBB). In stellar evolution models, the exact location of the RGBB depends on choices about the envelope mixing and undershoot parameters, which will likewise affect the structure and temperature at the base of the convection zone (Joyce & Chaboyer 2015;Khan et al. 2018).
Observational work has noted that several abundances, including lithium, [C/N], and 12 C/ 13 C change for stars near the red giant branch bump, particularly at low metallicity (see e.g. Carbon et al. 1982;Pilachowski 1986;Kraft 1994;Gratton et al. 2000). More recent work has made it clear that the amount of extra mixing in these stars is strongly metallicity dependent (Shetrone et al. 2019) and also depends on stellar mass (Magrini et al. 2021b).
Great theoretical efforts have been made to try to infer the physical cause of this extra mixing. Given that the extra mixing occurs in a region where the mean molecular weight gradient is inverted from 3 He burning (Eggleton et al. 2006), authors including Charbonnel & Zahn (2007b) pointed out that the region should be unstable to double diffusive instabilities, also known as thermohaline mixing. While there were challenges to this theory (Denissenkov 2010;Denissenkov & Merryfield 2011;Traxler et al. 2011;Brown et al. 2013;Sengupta & Garaud 2018;Garaud et al. 2019) on the grounds that the amount of mixing needed not did not quantitatively match the expectations given the fluid parameters, most authors have worked under the assumption that the extra mixing of all elements on the upper giant branch are somehow related to thermohaline mixing (e.g. Kirby et al. 2016;Charbonnel et al. 2020;Magrini et al. 2021a).
Previous theoretical explorations have indicated that the amount of extra mixing in this regime should depend on stellar mass and metallicity (e.g. Charbonnel & Lagarde 2010;Lagarde et al. 2017) and that the relevant process should affect several mixing tracers, including 12 C/ 13 C, [C/Fe], [N/Fe], and 7 Li, at approximately the same time above the red giant branch bump (Stancliffe et al. 2009;Charbonnel & Lagarde 2010). Predicting the exact amount of mixing or the resulting abundances is extremely challenging, however, as the process can be affected by rotation (Lagarde et al. 2011), magnetism (Charbonnel & Zahn 2007a), dynamical shear (Cantiello & Langer 2010), numerical model choices (Lattanzio et al. 2015), and other complications, as well as the detailed interactions, rather than just additions, of these processes (Maeder et al. 2013;Sengupta & Garaud 2018). Further, the estimated degree of im-pact of any one of these processes likewise depends on the stellar properties, as has been discussed extensively (see e.g. Cantiello & Langer 2010;Lagarde et al. 2011). However, by making reasonable physical assumptions, several authors have been able to construct predictions for the mixing that at least roughly match the observations (e.g. Lagarde et al. 2019;Henkel et al. 2017Henkel et al. , 2018Magrini et al. 2021b).
In the following manuscript, we present observational and theoretical evidence for distinct mixing causes for the depletion of lithium and [C/N], from which it follows that thermohaline mixing is not sufficient explanation for all of the extra mixing that occurs near the red giant branch bump. We first present simulations showing that the connection of the thermohaline mixing region with the surface convection zone corresponds to simultaneous depletion of both lithium and [C/N]. We follow with the presentation of observational data demonstrating that the abundance decreases in Li and [C/N] near the red giant branch bump are temporally distinct, with the depletion of lithium occurring earlier in a star's evolution than the depletion of [C/N] by more than 0.1 dex. We thus demonstrate that a theoretical framework ascribing these events to thermohaline mixing alone is inconsistent with the observations.

THEORETICAL EXPLORATION AND SIMULATIONS
Using the Modules for Experiments in Stellar Astrophysics (MESA; stable release version 15140) stellar evolution program, we evolve high structural and temporal resolution simulations of the interior mixing profiles of model stars in the approximate mass (0.8 − 1.2M ) and metallicity regime (−2 ≤ [Fe/H] ≤ 0 dex, or 0.0007 ≤ Z ≤ 0.02) of the observed sample (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019. We use the solar mixture of Grevesse & Sauval (1998) and OPAL opacities, metallicities corresponding to solar (Z = 0.02) and [Fe/H]= −1.2 (Z = 0.0009), the mixing length theory (MLT) implementation of Cox & Giuli (1968) with α MLT = 1.6, the Eddington T-τ grey model atmosphere for surface boundary conditions, the pp_extras.net nuclear reaction network, and the thermohaline mixing prescription of Kippenhahn et al. (1980) based on the analysis of Ulrich (1972) with a thermohaline efficiency coefficient α th = 2 (Eddington 1916;Böhm-Vitense 1958;Cox 1980;Iglesias & Rogers 1996a).The choice to set α th near unity was motivated by recent 2D and 3D hydrodynamic simulations showing that much lower values (e.g. α th = 1 vs α th = 667) are appropriate for stellar conditions as compared to salt fingers (Denissenkov 2010; Traxler et al. 2011;Brown et al. 2013). A value of α th = 2 specifically was chosen for consistency with Section 4.2 of MESA Instrument Paper II, in which MESA's thermohaline capabilities are described in more detail, and the associated test_suite case (Paxton et al. 2013). Using models described in Fraser et al. (2022), we find that variations of α th between 0.1 and 700 have no effect on our conclusions. The Ledoux criterion for convective stability is used throughout (Ledoux 1951;Anders et al. 2022). We improve MESA's default structural and temporal resolutions (mesh_delta_coeff and time_delta_coeff) by a factor of 10 to 100 in our calculations 1 . In physical units, our time steps in this regime correspond to 21,000 years-reduced to as low as 3000 years in convergence investigations-and our spatial resolution corresponds to roughly 5500 zones. These are above the thresholds necessary to resolve the relevant changes due to thermohaline mixing according to, for example, Lattanzio et al. (2015). We note, however, that while even higher resolutions might be necessary to argue for particular values of calibrated parameters (for example, the mixing efficiency); this is not necessary in our case, and we propose no such values. It is sufficient for our argument to demonstrate that neither time step nor structural resolution affect the direction of propagation of the thermohaline front or the simultaneity of dilution events. This is indeed the case (see, e.g. Fraser et al. 2022, which includes extensive resolution tests using an identical modeling configuration).
The models are terminated at roughly log g ≤ 1 on the upper red giant branch, after the RGBB. Figure 1 shows an evolutionary track from one such simulation, centered on the evolutionary regime where thermohaline mixing is hypothesized to be important. Locations in the HR diagram and values of log g are indicated on the track for six time steps dt corresponding to the structural snapshots presented subsequently. These are indexed beginning with one for convenience; in reality, there are thousands.
Our simulations trace the evolution of interior lithium, carbon, and nitrogen abundances and the nuclear energy generation front, ε nuc , as well as the locations and extents of various mixing regions. Figure 2 shows six snapshots of the interior near the base of the surface convection zone as a function of mass coordinate, m/M , with corresponding evolutionary locations indicated by equivalent dt on the HR diagram in Figure 1. The panels depict the surface-ward propagation and expansion of the thermohaline mixing front, shown in green, with 1 Full details of our configuration will be made publicly available on Zenodo upon publication. A basic template for simulating thermohaline mixing is freely available in the MESA test_suite.  Figure 1. The evolutionary track for our canonical thermohaline propagation simulation (1M , solar metallicity Z = 0.020) is shown, zoomed in on the RGBB. Surface gravity is indicated at six time steps, denoted dt. These index the evolutionary times at which the subsequent stellar interior snapshots are taken. The time step indices are revised to begin at one for the purpose of demonstration only.
increasing time: the region broadens and moves outward in mass coordinate from dt = 1 to dt = 4. At dt = 5, the thermohaline region connects with the surface convection zone, shown in pink, which is receding towards the surface at a much slower rate over this time interval (the recession of the surface convection zone is apparent at enhanced scale, shown in the following figure). At dt = 6, the regions have begun rapid mixing. Though panels dt = 5 and dt = 6 appear identical here, they represent summaries of the configuration before and after rapid mixing, respectively, as elaborated in Figure  3.
The black curve represents the nuclear burning shell, which also moves towards the surface over time. The blue and green curves show the interior abundance profiles of (C 12 /N 14 ) and lithium, respectively, where the lithium abundance has been shifted uniformly so as to appear on the same vertical scale as the other quantities. The local minimum of the lithium abundance near m/M = 0.25 in panel dt = 1 indicates the point at which the interior temperature is hot enough to destroy lithium. For detailed discussion on the behavior of the mean molecular weight profile and its inversion during the conditions that generate thermohaline mixing, see Charbonnel & Zahn (2007b) and Charbonnel & Lagarde (2010). Figure 3 shows the same six structural snapshots, but rescaled so that the changes in abundance are visually detectable. As the panels show, there is a slow decrease in [C 12 /N 14 ] and lithium as the time step increases from dt = 1 to dt = 5, but a sudden, much larger decrease in both at dt = 6. In the final panel, we see that the drops in lithium and [C/N] occur at the same instant to within model precision (taken to be an uncertainty of half of one time step). This is the instant at which the thermohaline front meets the surface convection zone. We also observe this directly in Figure 4, where we show the predicted abundances as a function of the surface gravity. In our models, the lithium and carbon-to-nitrogen ratios experience the onset of dilution at the same time step. The drop beginning at that point is significant relative to other changes in the abundances and should be visible in observational data. We argue that this conclusion is not sensitive to the numerical or physical choices made in the model. We emphasize, however, that absolute abundances and timings are extremely challenging to predict accurately given the wide variety of physical uncertainties in the models, (see e.g. Lattanzio et al. 2015, Cinquegrana et al., subm., for more extensive discussion). We are interested only in the relative changes here.
In general, this sequence is consistent with previous explorations of thermohaline mixing in low-mass red giant branch stars. After the RGBB, a region near the hydrogen burning shell becomes unstable to thermohaline mixing. On its own, thermohaline mixing in the model is insufficient to encourage mixing in the entire region between the burning shell and the base of the convection zone. However, when both rotational mixing and thermohaline mixing are included, the unstable region expands outwards, mixing through the region where lithium is destroyed, until it reaches the base of the surface convection zone. At this point, the star is able to (1) mix nitrogen enriched material from the shell up into the surface convection zone, and (2) to mix fragile elements like lithium from the cool surface of the star down into a region hot enough to destroy them. It is at this time that the surface abundances of lithium and [C/N] should change.

DATA
The theoretical framework discussed above makes a clear prediction that the abundances of all elements affected by thermohaline mixing should begin to change at the same point. The recent explosion of available spectroscopic measurements of well characterized stars makes it possible to test this prediction empirically.
For this analysis, we take the location of the red giant branch bump as a function of metallicity and the location of the change in the surface [C/N] ratio from  (Gunn et al. 2006) for a variety of reasons (Zasowski et al. 2017), but processed and analyzed homogenously by the ASPCAP pipeline (Nidever et al. 2015;Zamora et al. 2015;García Pérez et al. 2016;Holtzman et al. 2015Holtzman et al. , 2018. To separate red clump stars from the red giants used here, the authors used the spectroscopic evolutionary states estimated automatically from APOGEE (Holtzman et al. , 2018, which use a combination of temperature, gravity, and luminosity information, and are calibrated based on the sample of APOGEE stars with asteroseismic estimates of evolutionary state (Elsworth et al. 2017(Elsworth et al. , 2019. Since α-enhanced stars were formed when the galaxy was young, these stars are assumed to be old, and therefore consistently low mass (∼ 0.9 M ). Shetrone et al. (2019) therefore divided the full sample into metallicity bins, and studied each metallicity bin as if it represented an evolutionary sequence, roughly equivalent to a cluster. At each metallicity, they were able to identify the location of the overdensity of first ascent red giant branch stars that represented the red giant branch bump. In each bin, where possible, they also identified the location of the first dredge up, the extra mixing taking place near the red giant branch bump, and the location of any additional mixing on the upper giant branch if such mixing was detected at high significance (See Figure 5, or the original figures in that text for more detail). We show in Figure 6 their reported locations of the red giant branch bump and the extra mixing event, as probed by [C/N]. While the exact surface gravity of the red giant branch bump is metallicity dependent, the extra mixing event seems to happen consistently a few tenths of a dex above (after) the red giant branch bump. The decline of [C/N] at this point is consistent with the expectations from our theoretical understanding of thermohaline mixing as presented above.
We supplement these data with measurements of the lithium abundance of thousands of red giants in globular clusters observed by Kirby et al. (2016). From their dataset, we extract stars identified as first ascent red giants based on the inspection of their colors and magnitudes relative to other stars in the cluster (see Kirby et al. 2016) for more discussion.   Figure 3. Same as Figure 2 but zooming in on the abundances to emphasize that the abundances of lithium (green dot-dashed line) and (C/N) (purple solid line) change only slightly until the thermohaline region (green) meets the base of the surface convection zone (pink), at which point the abundances of both elements at the surface begin to drop together. The abundances are rescaled so as to appear on similar axes by linearly shifting the log 10 abundance of lithium by a factor of roughly 15. We emphasize that we do not intend or attempt to provide predictions for specific values of abundances in this study: that is a much more difficult problem, impacted by numerous physical assumptions and their uncertainties that are not relevant to our investigation of the sequence of dilution events alone.  The red vertical line marks the time step at which the dilution event related to the thermohaline front meeting the surface convection zone begins. This occurs slightly after the red giant branch bump, indicated with a purple dashed line. We show that this drop is predicted to begin at the same time step in both elements (emphasized in the zoomed-in panels on the left), and that the decline in the surface abundances at this point is significant compared to the slight changes visible earlier than and near the red giant branch bump (emphasized in the zoomed-out right panels). perbolic tangent function to the data between 1.5 and 2.9 dex in each metallicity bin . We report our results in Table 1. Unfortunately, there are insufficient data points  Nataf et al. (2013), and the V magnitude of their red giant branch bumps were measured. For those clusters, we used stars near that V magnitude with measured surface gravities to estimate the surface gravity of the red giant branch bump on the Kirby et al.
(2016) scale, and used the average metallicity measured by Kirby et al. (2016) for the whole cluster. While there may be uncertainties and systematics from this method for individual clusters, the ensemble seems to match well the trend in bump surface gravity measured in Shetrone et al. (2019) (Figure 6). We therefore make the assumption that the two measurements are on the same scale and can be compared, a point to which we return in the discussion section. We show in Figure 6 that a decrease in lithium abundance does indeed happen near the red giant branch bump, but it is not coincident with the drop in [C/N]. This implies that the two changes cannot be due to the same thermohaline mixing event.

DISCUSSION
In this paper, we have presented an argument regarding the the observed timing of chemical depletion events and the ability of thermohaline mixing alone to explain them. Our simulations clearly show it cannot, thus inviting a more sophisticated theoretical explanation for the observed time-separation of abundance depletion trends.
Given the observational data available, the extra mixing events that alter the surface abundances of [C/N] and lithium cannot both be ascribed to thermohaline mixing since they do not occur at the same time. The regions that must be mixed have previously been homogenized during the first dredge up, and therefore should not retain chemical gradients and discontinuities left behind from pre-main-sequence or main sequence processes. As shown in the literature and by the modeling above, the thermohaline mixing front is propagating outwards from the core after the bump, and therefore extra mixing events related to thermohaline mixing should not begin to alter the surface chemistry until about 0.1 dex above the red giant branch bump. While this result is consistent with the measured location of the change in [C/N], it is inconsistent with the location of the change of the lithium abundance. Such discrepancies have been noted in the past in detailed cluster analyses (Angelou et al. 2015), but in general it has been argued that they result from other other issues in the analysis such as the choice of distance modulus (Henkel et al. 2017). Under the assumption that the offset we see in our data is real, we first discuss possible theoretical explanations for this result, and then suggest several observational tests which could ensure that it is robust.
Theoretically, in order to obtain a drop in lithium before the thermohaline front has met the base of the convection zone, there needs to be sufficient extra mixing beneath the convective envelope to penetrate into the region of the star that is hot enough to destroy lithium. Simulations of convective zones do often see convective-like flows moving into regions that are formally stably stratified, variously called convective overshoot (e.g. Zahn 1991; Korre et al. 2019), entrainment (e.g. Meakin & Arnett 2007;Horst et al. 2021) or convective penetration (e.g. Hurlburt et al. 1994;Anders et al. 2021). However, both simulations and observations (Claret & Torres 2019;Pedersen et al. 2021) of main sequence and red giant branch stars generally find these overshoot regions to be significantly smaller than a pressure scale height (Joyce & Chaboyer 2018). In the case of red giant stars, the lithium burning zone is more than a pressure scale height below the formal boundary of the surface convection zone (Straus et al. 1976). Rotational mixing can be present in this region in models (see e.g. Lagarde et al. 2012;Charbonnel et al. 2020), but it is generally not expected to be strong enough to deplete a significant amount of lithium on its own. This suggests that some unusual process, or combination of processes, may be acting in this region.
However, we also acknowledge the possibility that the offset between the lithium and [C/N] mixing locations could be the result of observational biases. In particular, calibrating the absolute surface gravities and metallicities between surveys can be challenging, and so exploration of a sample of stars where both lithium and [C/N] have been measured would be a strong test of this result, as would detailed analysis of both elements in a single cluster, where systematics are easier to control and differential measurements are extremely powerful. Larger sample sizes, from, e.g., the more recent APOGEE DR17 data or GALAH DR3, would also reduce the observational uncertainty in the location of the drop in abundance for each element, as well as the location and size of the luminosity bump. More extensive catalogs might also be better able to map the shape of the change near the dilution drop. They may also be better able to account for the expected differences in the rate of change of lithium versus the carbon-to-nitrogen ratio, which could affect the inferred location of the onset of the mixing. Additional mixing diagnostics such as 12 C/ 13 C might also be illuminating, as would a study of the mass dependence of the location and depth of extra mixing. Mis-identification of evolutionary states in the data could also alter the measured location of the extra mixing: for example, dragging the estimated drop from the bump towards the clump, where stars are already mixed. This would present as a shift towards higher gravity at low metallicity and towards a lower gravity at high metallicity. While we believe that this is not a significant factor for our data sets, where evolutionary states have been carefully calibrated or curated by eye, it is important to consider for future investigations of this type.
This investigation serves as an exciting foray into the new reality where the properties of extremely large samples of stars are exquisitely characterized, and therefore even minute surface changes can be identified observationally and used to test theories of stellar interiors and evolution. As spectroscopic surveys release even more data (e.g. GALAH, Gaia-ESO, SDSS-V), and as these data can be combined with other methods of characterizing stars, including astrometry, asteroseismology, interferometry, photometry, and time domain studies, we venture into a regime where a theory must not only qualitatively match the expected results, but also quantitatively agree with precisely characterized observational data, with significantly reduced systematic uncertainties.