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

Jamie Tayar These authors contributed equally to this work  Department of
Astronomy, University of Florida, Bryant Space Science Center, Stadium Road,
Gainesville, FL 32611, USA  Institute for Astronomy, University of Hawai‘i at
Mānoa, 2680 Woodlawn Drive, Honolulu, HI 96822, USA  NASA Hubble Fellow
Meridith Joyce These authors contributed equally to this work  Space Telescope
Science Institute, 3700 San Martin Dr, Baltimore, MD 21202, USA  Kavli
Institute for Theoretical Physics, University of Santa Barbara, California,
93106, USA Jamie Tayar<EMAIL_ADDRESS>Meridith Joyce<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

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.

stellar evolution, stellar abundances, stellar interiors, thermohaline mixing

##  1 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 towards 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
12C/13C 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 structure
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 12C/13C 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 3He 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 12C/13C, [C/Fe],
[N/Fe], and 7Li, 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 impact 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., 2017, 2018; Magrini 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.

##  2 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_{\odot}$) and metallicity regime
($-2\leq$ [Fe/H] $\leq 0$ dex, or $0.0007\leq Z\leq$ 0.02) of the observed
sample (Paxton et al., 2011, 2013, 2015, 2018, 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 $\alpha_{\text{MLT}}=1.6$, the
Eddington T-$\tau$ 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 $\alpha_{\text{th}}=2$
(Eddington, 1916; Böhm-Vitense, 1958; Cox, 1980; Iglesias & Rogers, 1996a).The
choice to set $\alpha_{\text{th}}$ near unity was motivated by recent 2D and
3D hydrodynamic simulations showing that much lower values (e.g.
$\alpha_{\text{th}}=1$ vs $\alpha_{\text{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 $\alpha_{\text{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 $\alpha_{\text{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
calculations111Full 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`.. 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\leq 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.

Figure 1: The evolutionary track for our canonical thermohaline propagation
simulation ($1M_{\odot}$, 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.

Our simulations trace the evolution of interior lithium, carbon, and nitrogen
abundances and the nuclear energy generation front,
$\varepsilon_{\text{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_{\star}$,
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 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 (C12/N14) 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_{\star}=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 [C12/N14] 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.

Figure 2: Structural snapshots at evolutionary time steps for the
$1M_{\odot}$, solar composition model, as indicated in Figure 1, are shown. As
the model steps through time (marked in the upper right corner) the nuclear
burning shell (black dashed line) moves outwards, as the surface convection
zone’s (pink shaded region) size is reduced. The thermohaline region (green
shaded region) grows from the hydrogen burning shell until it meets the base
of the convective region. The abundance of lithium (green line) and the ratio
of carbon to nitrogen (blue line) are also shown, with units rescaled to fit
within the plot. In this simulation, the lithium depletion region and the
(C/N) depleted region are connected to the surface convection zone by the
thermohaline front at the same time.

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.

Figure 4: Model predictions for the abundances of lithium (top) and [C/N]
(bottom) are shown as the star evolves in surface gravity. 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).

##  3 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 Shetrone et al. (2019). The authors used a large dataset of
$\alpha$-element enhanced stars from the SDSS-IV’s (Blanton et al., 2017)
APOGEE-2 (Majewski et al., 2017) Data Release 14 (Abolfathi et al., 2018).
These stars were observed by the APOGEE (Wilson et al., 2019) spectrograph on
the Sloan Foundation Telescope (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., 2015, 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., 2015, 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, 2019).

Figure 5: With large samples of stars binned in metallicity, it is possible to
trace chemical evolution as a function of surface gravity using lithium (teal
diamonds) and [C/N] (blue triangles) as well as empirically identify the red
giant branch bump (purple circle, and descending line for reference). In
analogy with the figures in Shetrone et al. (2019) and in order to make
changes in abundance more clear, we have indicated with lines the approximate
level of the post-dredge-up, pre-bump abundances (teal and blue lines).

Since $\alpha$-enhanced stars were formed when the galaxy was young, these
stars are assumed to be old, and therefore consistently low mass ($\sim 0.9$
$\mathrm{M_{\odot}}$). 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.

Table 1: Locations on the giant branch for the drop in surface lithium abundance measured from the Kirby et al. (2016) data (see Section 3) as well as the drop in [C/N] and the red giant branch bump as measured in Shetrone et al. (2019). Error bars represent the uncertainty in the measured location from various methods rather than the systematic uncertainties on the data. [M/H] | log(g)Li | log(g)CN | log(g)bump
---|---|---|---
0.0 | … | … | 2.54$\pm$0.01
-0.2 | … | … | 2.56$\pm$0.01
-0.4 | … | 2.28$\pm$0.08 | 2.55$\pm$0.04
-0.6 | … | 2.32$\pm$0.11 | 2.52$\pm$0.12
-0.8 | … | 2.36$\pm$0.11 | 2.46$\pm$0.01
-1.0 | 2.37$\pm$0.15 | 2.23$\pm$0.09 | 2.30$\pm$0.04
-1.2 | 2.42$\pm$0.19 | 2.14$\pm$0.03 | 2.33$\pm$0.05
-1.4 | 2.22$\pm$0.21 | 2.14 $\pm$0.03 | 2.25$\pm$0.01
-1.6 | 2.15$\pm$0.15 | … | 2.27$\pm$0.07
-1.8 | 2.47$\pm$0.10 | … | …
-2.0 | 2.36$\pm$0.06 | … | …
-2.2 | 2.28$\pm$0.01 | … | …
-2.4 | 1.98$\pm$0.01 | … | …
-2.6 | 2.22$\pm$0.22 | … | …

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. Following the
procedure described in Shetrone et al. (2019), we bin the stars in metallicity
and search for a drop in lithium near the red giant branch bump (see data in
Figure 5) by fitting a hyperbolic 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 in this sample to robustly
detect the red giant branch bump location directly in the Kirby et al. (2016)
spectroscopic sample. However some, of these clusters were studied with HST by
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.

Figure 6: The location in surface gravity, as a function of metallicity, of
the drop in [C/N] associated with extra mixing (blue triangles), and the drop
in lithium near the same point (teal diamonds). The empirical locations of the
red giant branch bump from Shetrone et al. (2019) are marked in filled
circles, and the estimated locations of the bump in Kirby et al. (2016)
clusters that were studied by Nataf et al. (2013) are marked as open circles.
While the drop in [C/N] is always above the red giant branch bump as expected,
the drop in lithium abundance seems to happen significantly earlier,
suggesting that a different mechanism is responsible.

##  4 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 12C/13C 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.

## Acknowledgements

J.T. and M.J. contributed equally to this manuscript. The authors gratefully
acknowledge helpful discussions with Adrian Fraser, Matteo Cantiello, Evan
Anders, and the rest of the thermohaline working group at the KITP Probes of
Transport in Stars 2021 meeting. J.T. would like to thank Evan Kirby for
discussions that contributed to the beginnings of this work and Ann Boesgaard
for pointing out useful references. M.J. wishes to thank Brian Chaboyer for
useful discussion on numerics during the second revision of this manuscript.
J.T. and M.J. also wish to thank Amanda Karakas, Marc Pinsonneault, and László
Molnár, each of whom provided an independent reading and review of this
manuscript. We thank the referees for their comments that improved this
manuscript.

J.T. acknowledges that support for this work was provided by NASA through the
NASA Hubble Fellowship grant No.51424 awarded by the Space Telescope Science
Institute, which is operated by the Association of Universities for Research
in Astronomy, Inc., for NASA, under contract NAS5-26555. M.J. was supported by
the Lasker Data Science Fellowship, awarded by the Space Telescope Science
Institute, and the Kavli Institute for Theoretical Physics at the University
of California, Santa Barbara. This research was supported in part by the
National Science Foundation under Grant No. NSF PHY-1748958.

##  Appendix A MESA Details

The MESA EOS is a blend of the OPAL (Rogers & Nayfonov, 2002), SCVH (Saumon et
al., 1995), FreeEOS (Irwin, 2004), HELM (Timmes & Swesty, 2000), and PC
(Potekhin & Chabrier, 2010) EOSes.

Radiative opacities are primarily from OPAL (Iglesias & Rogers, 1993, 1996b),
with low-temperature data from Ferguson et al. (2005) and the high-
temperature, Compton-scattering dominated regime by Buchler & Yueh (1976).
Electron conduction opacities are from Cassisi et al. (2007).

Nuclear reaction rates are from JINA REACLIB (Cyburt et al., 2010) plus
additional tabulated weak reaction rates Fuller et al. (1985); Oda et al.
(1994); Langanke & Martínez-Pinedo (2000). Screening is included via the
prescription of Chugunov et al. (2007). Thermal neutrino loss rates are from
Itoh et al. (1996).

## References

  * Abolfathi et al. (2018) Abolfathi, B., Aguado, D. S., Aguilar, G., et al. 2018, ApJS, 235, 42
  * Anders et al. (2021) Anders, E. H., Jermyn, A. S., Lecoanet, D., & Brown, B. P. 2021, arXiv e-prints, arXiv:2110.11356
  * Anders et al. (2022) Anders, E. H., Jermyn, A. S., Lecoanet, D., et al. 2022, arXiv e-prints, arXiv:2203.06186
  * Angelou et al. (2015) Angelou, G. C., D’Orazi, V., Constantino, T. N., et al. 2015, MNRAS, 450, 2423
  * Blanton et al. (2017) Blanton, M. R., Bershady, M. A., Abolfathi, B., et al. 2017, AJ, 154, 28
  * Boesgaard et al. (2020) Boesgaard, A. M., Lum, M. G., & Deliyannis, C. P. 2020, ApJ, 888, 28
  * Böhm-Vitense (1958) Böhm-Vitense, E. 1958, ZAp, 46, 108
  * Boothroyd & Sackmann (1999) Boothroyd, A. I., & Sackmann, I. J. 1999, ApJ, 510, 232
  * Brown et al. (2013) Brown, J. M., Garaud, P., & Stellmach, S. 2013, ApJ, 768, 34
  * Buchler & Yueh (1976) Buchler, J. R., & Yueh, W. R. 1976, ApJ, 210, 440
  * Cantiello & Langer (2010) Cantiello, M., & Langer, N. 2010, A&A, 521, A9
  * Carbon et al. (1982) Carbon, D. F., Romanishin, W., Langer, G. E., et al. 1982, ApJS, 49, 207
  * Cassisi et al. (2007) Cassisi, S., Potekhin, A. Y., Pietrinferni, A., Catelan, M., & Salaris, M. 2007, ApJ, 661, 1094
  * Chanamé et al. (2005) Chanamé, J., Pinsonneault, M., & Terndrup, D. M. 2005, ApJ, 631, 540
  * Charbonnel & Lagarde (2010) Charbonnel, C., & Lagarde, N. 2010, A&A, 522, A10
  * Charbonnel et al. (1994) Charbonnel, C., Vauclair, S., Maeder, A., Meynet, G., & Schaller, G. 1994, A&A, 283, 155
  * Charbonnel & Zahn (2007a) Charbonnel, C., & Zahn, J. P. 2007a, A&A, 476, L29
  * Charbonnel & Zahn (2007b) —. 2007b, A&A, 467, L15
  * Charbonnel et al. (2020) Charbonnel, C., Lagarde, N., Jasniewicz, G., et al. 2020, A&A, 633, A34
  * Chugunov et al. (2007) Chugunov, A. I., Dewitt, H. E., & Yakovlev, D. G. 2007, Phys. Rev. D, 76, 025028
  * Claret & Torres (2019) Claret, A., & Torres, G. 2019, ApJ, 876, 134
  * Cox (1980) Cox, J. P. 1980, Theory of stellar pulsation
  * Cox & Giuli (1968) Cox, J. P., & Giuli, R. T. 1968, Principles of stellar structure
  * Cyburt et al. (2010) Cyburt, R. H., Amthor, A. M., Ferguson, R., et al. 2010, ApJS, 189, 240
  * Denissenkov (2010) Denissenkov, P. A. 2010, ApJ, 723, 563
  * Denissenkov & Merryfield (2011) Denissenkov, P. A., & Merryfield, W. J. 2011, ApJ, 727, L8
  * Eddington (1916) Eddington, A. S. 1916, MNRAS, 77, 16
  * Eggleton et al. (2006) Eggleton, P. P., Dearborn, D. S. P., & Lattanzio, J. C. 2006, Science, 314, 1580
  * Elsworth et al. (2017) Elsworth, Y., Hekker, S., Basu, S., & Davies, G. R. 2017, MNRAS, 466, 3344
  * Elsworth et al. (2019) Elsworth, Y., Hekker, S., Johnson, J. A., et al. 2019, MNRAS, 489, 4641
  * Ferguson et al. (2005) Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585
  * Fraser et al. (2022) Fraser, A. E., Joyce, M., Anders, E. H., Tayar, J., & Cantiello, M. 2022, arXiv e-prints, arXiv:2204.08487
  * Fuller et al. (1985) Fuller, G. M., Fowler, W. A., & Newman, M. J. 1985, ApJ, 293, 1
  * Fusi Pecci et al. (1990) Fusi Pecci, F., Ferraro, F. R., Crocker, D. A., Rood, R. T., & Buonanno, R. 1990, A&A, 238, 95
  * Garaud et al. (2019) Garaud, P., Kumar, A., & Sridhar, J. 2019, ApJ, 879, 60
  * García Pérez et al. (2016) García Pérez, A. E., Allende Prieto, C., Holtzman, J. A., et al. 2016, AJ, 151, 144
  * Gratton et al. (2000) Gratton, R. G., Sneden, C., Carretta, E., & Bragaglia, A. 2000, A&A, 354, 169
  * Grevesse & Sauval (1998) Grevesse, N., & Sauval, A. J. 1998, Space Sci. Rev., 85, 161
  * Gunn et al. (2006) Gunn, J. E., Siegmund, W. A., Mannery, E. J., et al. 2006, AJ, 131, 2332
  * Henkel et al. (2018) Henkel, K., Karakas, A. I., Casey, A. R., Church, R. P., & Lattanzio, J. C. 2018, ApJ, 863, L5
  * Henkel et al. (2017) Henkel, K., Karakas, A. I., & Lattanzio, J. C. 2017, MNRAS, 469, 4600
  * Holtzman et al. (2015) Holtzman, J. A., Shetrone, M., Johnson, J. A., et al. 2015, AJ, 150, 148
  * Holtzman et al. (2018) Holtzman, J. A., Hasselquist, S., Shetrone, M., et al. 2018, AJ, 156, 125
  * Horst et al. (2021) Horst, L., Hirschi, R., Edelmann, P. V. F., Andrássy, R., & Röpke, F. K. 2021, A&A, 653, A55
  * Hurlburt et al. (1994) Hurlburt, N. E., Toomre, J., Massaguer, J. M., & Zahn, J.-P. 1994, ApJ, 421, 245
  * Iben (1964) Iben, Jr., I. 1964, ApJ, 140, 1631
  * Iben (1967) —. 1967, ApJ, 147, 624
  * Iglesias & Rogers (1993) Iglesias, C. A., & Rogers, F. J. 1993, ApJ, 412, 752
  * Iglesias & Rogers (1996a) —. 1996a, ApJ, 464, 943
  * Iglesias & Rogers (1996b) —. 1996b, ApJ, 464, 943
  * Irwin (2004) Irwin, A. W. 2004, The FreeEOS Code for Calculating the Equation of State for Stellar Interiors I: An Improved EFF-Style Approximation for the Fermi-Dirac Integrals
  * Itoh et al. (1996) Itoh, N., Hayashi, H., Nishikawa, A., & Kohyama, Y. 1996, ApJS, 102, 411
  * Joyce & Chaboyer (2015) Joyce, M., & Chaboyer, B. 2015, ApJ, 814, 142
  * Joyce & Chaboyer (2018) —. 2018, ApJ, 856, 10
  * Khan et al. (2018) Khan, S., Hall, O. J., Miglio, A., et al. 2018, ApJ, 859, 156
  * Kippenhahn et al. (1980) Kippenhahn, R., Ruschenplatt, G., & Thomas, H. C. 1980, A&A, 91, 175
  * Kippenhahn & Weigert (1994) Kippenhahn, R., & Weigert, A. 1994, Stellar Structure and Evolution
  * Kirby et al. (2016) Kirby, E. N., Guhathakurta, P., Zhang, A. J., et al. 2016, ApJ, 819, 135
  * Korre et al. (2019) Korre, L., Garaud, P., & Brummell, N. H. 2019, MNRAS, 484, 1220
  * Kraft (1994) Kraft, R. P. 1994, PASP, 106, 553
  * Lagarde et al. (2011) Lagarde, N., Charbonnel, C., Decressin, T., & Hagelberg, J. 2011, A&A, 536, A28
  * Lagarde et al. (2012) Lagarde, N., Decressin, T., Charbonnel, C., et al. 2012, A&A, 543, A108
  * Lagarde et al. (2017) Lagarde, N., Robin, A. C., Reylé, C., & Nasello, G. 2017, A&A, 601, A27
  * Lagarde et al. (2019) Lagarde, N., Reylé, C., Robin, A. C., et al. 2019, A&A, 621, A24
  * Lambert (1981) Lambert, D. L. 1981, in Astrophysics and Space Science Library, Vol. 88, Physical Processes in Red Giants, ed. J. Iben, I. & A. Renzini, 115–134
  * Langanke & Martínez-Pinedo (2000) Langanke, K., & Martínez-Pinedo, G. 2000, Nuclear Physics A, 673, 481
  * Lattanzio et al. (2015) Lattanzio, J. C., Siess, L., Church, R. P., et al. 2015, MNRAS, 446, 2673
  * Ledoux (1951) Ledoux, P. 1951, ApJ, 114, 373
  * Maeder et al. (2013) Maeder, A., Meynet, G., Lagarde, N., & Charbonnel, C. 2013, A&A, 553, A1
  * Magrini et al. (2021a) Magrini, L., Smiljanic, R., Lagarde, N., et al. 2021a, The Messenger, 185, 18
  * Magrini et al. (2021b) Magrini, L., Lagarde, N., Charbonnel, C., et al. 2021b, A&A, 651, A84
  * Majewski et al. (2017) Majewski, S. R., Schiavon, R. P., Frinchaboy, P. M., et al. 2017, AJ, 154, 94
  * Martig et al. (2016) Martig, M., Fouesneau, M., Rix, H.-W., et al. 2016, MNRAS, 456, 3655
  * Masseron & Gilmore (2015) Masseron, T., & Gilmore, G. 2015, MNRAS, 453, 1855
  * Meakin & Arnett (2007) Meakin, C. A., & Arnett, D. 2007, ApJ, 667, 448
  * Nataf et al. (2013) Nataf, D. M., Gould, A. P., Pinsonneault, M. H., & Udalski, A. 2013, ApJ, 766, 77
  * Ness et al. (2016) Ness, M., Hogg, D. W., Rix, H.-W., et al. 2016, ApJ, 823, 114
  * Nidever et al. (2015) Nidever, D. L., Holtzman, J. A., Allende Prieto, C., et al. 2015, AJ, 150, 173
  * Oda et al. (1994) Oda, T., Hino, M., Muto, K., Takahara, M., & Sato, K. 1994, Atomic Data and Nuclear Data Tables, 56, 231
  * Palacios et al. (2006) Palacios, A., Charbonnel, C., Talon, S., & Siess, L. 2006, A&A, 453, 261
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10
  * Pedersen et al. (2021) Pedersen, M. G., Aerts, C., Pápics, P. I., et al. 2021, Nature Astronomy, arXiv:2105.04533
  * Pilachowski (1986) Pilachowski, C. 1986, ApJ, 300, 289
  * Potekhin & Chabrier (2010) Potekhin, A. Y., & Chabrier, G. 2010, Contributions to Plasma Physics, 50, 82
  * Rogers & Nayfonov (2002) Rogers, F. J., & Nayfonov, A. 2002, ApJ, 576, 1064
  * Saumon et al. (1995) Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713
  * Sengupta & Garaud (2018) Sengupta, S., & Garaud, P. 2018, ApJ, 862, 136
  * Shetrone et al. (2019) Shetrone, M., Tayar, J., Johnson, J. A., et al. 2019, ApJ, 872, 137
  * Stancliffe et al. (2009) Stancliffe, R. J., Church, R. P., Angelou, G. C., & Lattanzio, J. C. 2009, MNRAS, 396, 2313
  * Straus et al. (1976) Straus, J. M., Blake, J. B., & Schramm, D. N. 1976, ApJ, 204, 481
  * Sweigart et al. (1989) Sweigart, A. V., Greggio, L., & Renzini, A. 1989, ApJS, 69, 911
  * Timmes & Swesty (2000) Timmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501
  * Traxler et al. (2011) Traxler, A., Garaud, P., & Stellmach, S. 2011, ApJ, 728, L29
  * Ulrich (1972) Ulrich, R. K. 1972, ApJ, 172, 165
  * Wilson et al. (2019) Wilson, J. C., Hearty, F. R., Skrutskie, M. F., et al. 2019, PASP, 131, 055001
  * Zahn (1991) Zahn, J. P. 1991, A&A, 252, 179
  * Zamora et al. (2015) Zamora, O., García-Hernández, D. A., Allende Prieto, C., et al. 2015, AJ, 149, 181
  * Zasowski et al. (2017) Zasowski, G., Cohen, R. E., Chojnowski, S. D., et al. 2017, AJ, 154, 198