Observations of globular clusters and low-metallicity field stars show significant changes in the abundance ratios of elements known to be sensitive to mixing, including 12 C/ 13 C, lithium, and [C/N], as a star evolves up the red giant branch (Carbon et al. 1982;Pilachowski 1986;Kraft 1994;Shetrone et al. 2019). These changes occur around the red giant branch bump (RGBB) and are largest in the most metal-poor stars (e.g. Gratton et al. 2000). Large samples of stars that can be used to trace this mixing are now available from a variety of spectroscopic surveys, including GALAH (Buder et al. 2019), APOGEE (Abdurro'uf et al. 2022), and GAIA-ESO (Magrini et al. 2021a). However, observed surface abundance trends are in tension with standard theoretical stellar evolution models, which predict that the surface chemistry should not evolve in this regime.

As low-mass stars ascend the red giant branch, they undergo a series of mixing and homogenizing events as their interior burning and energy transport zones interact. Near the base of the red giant branch, the surface convection zone reaches its deepest level of penetration into the stellar interior, leaving behind a chemical discontinuity from which it recedes in subsequent evolution. This inflection in the convection zone's movement is known as the "first dredge-up." The red giant branch bump occurs when the outward-propagating hydrogen burning shell encounters this chemical discontinuity, triggering a structural realignment in which the star's core contracts and the luminosity drops, causing a disruption to the otherwise monotonic increase in luminosity along the red giant branch (Christensen-Dalsgaard 2015). In one dimensional (1D) stellar models, the sensitivity of the RGBB to physical assumptions makes it a powerful diagnostic of interior mixing processes (e.g. Fusi Pecci et al. 1990;Salaris et al. 2002;Joyce & Chaboyer 2015;Khan et al. 2018). However, in standard models of red giant stars, there is no mixing between the hydrogen-burning shell and the overlying convective envelope after the first dredge-up, and no change in surface abundances is predicted in this regime. This is in direct conflict with abundance trends found in observations.

The most widely studied candidate mechanism for rectifying this discrepancy is thermohaline mixing, driven by an inversion of the mean molecular weight mu caused by 3 He burning. The potential for 3 He burning to cause a µ inversion and drive thermohaline mixing (though not specifically in low-mass RGB stars) was first pointed out by Ulrich (1972). Later, this µ inversion was identified as a significant source of mixing in low-mass RGB stars (Dearborn et al. 2006;Eggleton et al. 2006Eggleton et al. , 2007)), but was connected to the Rayleigh-Taylor instability, not thermohaline mixing. The specific connection between extra mixing in low-mass RGB stars and thermohaline mixing was later made by Charbonnel & Zahn (2007). As the hydrogen-burning shell moves into the region chemically homogenized by the first dredge-up, the 3 He( 3 He, 2p) 4 He reaction creates an inversion of the mean molecular weight µ. While this µ inversion is insufficient to generate a convective region (c.f. Cantiello & Langer 2010), these conditions give rise to the thermohaline instability, a phenomenon perhaps best known in the context of salt water in Earth's oceans (Stern 1960;Baines & Gill 1969).

Thermohaline mixing occurs in Ledoux-stable regions that have stably stratified temperature gradients but unstable mean molecular weight stratification (see, e.g., Garaud 2018 for a full review or Salaris & Cassisi 2017 for discussion focused on a 1D stellar modeling context). Thermohaline mixing is a double-diffusive phenomenon present in fluids that have different diffusivities for heat and chemical composition that, in turn, make opposing contributions to the radial density gradient (Turner 1974). This process may facilitate the radial mixing of elements between the hydrogen-burning shell and the stellar convective envelope, thus producing measurable changes in the surface mixing diagnostics after the first dredge-up.

Fluid dynamicists have studied thermohaline mixing in great detail, often employing 3D simulations in "local" domains, i.e., periodic boundary conditions with constant linear background gradients in temperature T and mean molecular weight µ under the Boussinesq approximation (Spiegel & Veronis 1960). Within this standard framework, the efficiency of thermohaline mixing, D th (the degree to which thermohaline mixing enhances chemical mixing via turbulent motions, as explained in Sec. 3 below), strictly depends on three dimensionless numbers, which we introduce here but describe in further detail in Secs. 2 and 3. The Prandtl number Pr and diffusivity ratio τ characterize molecular diffusivities and are fluid properties, meaning they are independent of background gradients. Pr is the ratio of the kinematic viscosity to thermal diffusivity, and τ is the ratio of the compositional and thermal diffusivities. In stars, Pr ≈ τ 1, so double-diffusive instabilities like thermohaline mixing are readily driven (Garaud 2018).

The third dimensionless quantity describing thermohaline mixing is the density ratio R 0 , which is the ratio of the temperature gradient's stabilizing contribution to the density divided by the µ gradient's destabilizing contribution to the density. The density ratio is a measure of how conducive to driving thermohaline mixing the the conditions in a radial location within a star are. If the density ratio is large, then the destabilizing µ gradient is weak and the system may be stable (if R 0 exceeds a threshold that depends on τ ) or only weakly unstable to the thermohaline instability. A small density ratio, on the other hand, corresponds to a significant inversion of the µ profile and thus the system is strongly unstable to thermohaline mixing and possibly even convection.

In this work, we study the reduced density ratio, r, which combines R 0 and τ into a single quantity that directly determines the instability of a system to the thermohaline instability (Traxler et al. 2011). The reduced density ratio directly defines the fluid dynamical stability of a system such that

System is convectively unstable

We stress that r is not a measure of the efficiency or mixing rate of thermohaline mixing. Rather, it is a measure of the structural stability of a system, or its tendency to mix. While mixing efficiency, D th , does depend on r, the exact dependence is an open question. Several simplified mixing prescriptions exist for predicting how efficiency depends on different fluid parameters (see review by Garaud 2018). As we will show in Sec. 3, often these mixing prescriptions diverge significantly in their predictions for how mixing efficiency (D th ) varies with r.

Given that the physical conditions required to trigger the thermohaline instability are in place at around the same time that extra mixing has been observed in red giant stars (e.g. Lagarde et al. 2015), most authors have assumed that all of the observed extra mixing can be attributed to the thermohaline instability (e.g. Kirby et al. 2016;Charbonnel et al. 2020;Magrini et al. 2021b). However, this connection has also been questioned for a number of reasons.

First, reproducing the observed amounts of mixing in this regime with 1D models requires assuming that thermohaline mixing is much more efficient than most fluid simulations would suggest is reasonable (Denissenkov 2010; Denissenkov & Merryfield 2011;Traxler et al. 2011;Brown et al. 2013). Questions have likewise been raised about whether the evolutionary timing of the observed extra mixing is truly consistent with thermohaline models (see e.g. Angelou et al. 2015;Henkel et al. 2017;Tayar & Joyce 2022).

Additionally, authors have put forth many different prescriptions for including thermohaline mixing in 1D models. Many previous works (for example, Charbonnel & Lagarde 2010, Lagarde et al. 2017) have used observations to calibrate or constrain 1D stellar model parameters that control the overall mixing efficiency within a particular mixing prescription (for example, α th in the Kippenhahn prescription, described in Sec. 3 below). These efforts generally employ a single prescription (e.g. that of Kippenhahn et al. 1980) and calibrate that prescription to observations. However, calibrating free parameters (e.g. α th , see Sec. 3) within any given mixing prescription, while useful in a 1D context, does not allow us to discriminate between the different proposed mixing prescriptions, which predict different (even, in some cases, opposite) trends in mixing efficiency D th versus fluid parameters like r. We are unaware of any instance in which observations of mixing have been used to (in)validate proposed theoretical mixing prescriptions.

This issue is all the more pressing in light of recent work demonstrating that models of thermohaline instability that include the presence of a relatively lowamplitude magnetic field can result in much more efficient mixing (larger D th ) and significantly different dependence of that efficiency on r (Harrington & Garaud 2019). Thus, while the issue that 1D models need to assume extremely efficient mixing in order to agree with observations could potentially be addressed by the inclusion of magnetic effects, it is not clear whether the profoundly different trend in mixing efficiency versus r is consistent with observations.

With more observational data available in stars across different masses and metallicities, a natural question to ask is: can these data be used to identify which of the different relationships between mixing efficiency (D th ) and r are more consistent with observations? This question is challenging because observations show how the amount of mixing mixing varies with mass and metallicity, while different mixing prescriptions relate mixing efficiency to r. Fluid properties (e.g. Pr, τ ) are readily extracted from 1D stellar evolution models (e.g. Jermyn et al. 2022), thus allowing values to be inferred for stars across different masses and metallicities. In contrast, we are unaware of any prior work which has extracted r from 1D stellar evolution models, which is required in order to understand how observed mixing trends relate to thermohaline models developed from 3D dynamical simulations.

Given these observational and theoretical questions, the development of a framework through which we can determine whether particular thermohaline prescriptions are more or less consistent with observations is timely and imperative. In this paper, we put forth such a framework: one that shifts the focus away from the calibration of individual mixing efficiency parameters within particular 1D thermohaline mixing prescriptions, and instead facilitates inter-comparison among the mixing models themselves.

We demonstrate a robust and modelagnostic means of relating the non-dimensional fluid parameters relevant to thermohaline mixing (particularly r) to the observed mixing around the RGB bump and show that this correlation is indeed qualitatively consistent with 1D prescriptions of thermohaline mixing informed by 3D simulations.

Further, while previous work (e.g. Charbonnel & Zahn 2007) has used the measurements of the overall amount of extra mixing to tune model parameters that control the overall efficiency of thermohaline mixing prescriptions, our framework allows us to use trends in extra mixing as a function of fundamental stellar parameters to probe trends in mixing efficiency (D th ) as a function of fluid parameters. Using this framework, we demonstrate that the magnetized thermohaline mixing prescription put forth by Harrington & Garaud (2019), which addressed significant outstanding problems by enhancing mixing efficiency over hydrodynamic levels, is not consistent with observations (when fixing their magnetic parameter H B ).

This paper is organized as follows: we begin by summarizing the formalism and stellar structure quantities relevant to thermohaline mixing (Sec. 2). This is followed by a description of various 1D mixing prescriptions commonly adopted in stellar evolution calculations (Sec. 3). We then introduce a suite of 1D MESA simulations and calculate the relevant fluid parameters in the thermohaline region for a range of mass and metallicity assumptions (Secs. 4 and 5). Finally, we compare an observational proxy of extra mixing, the decrease in [C/N] near the RGBB, to theoretical trends predicted by existing 1D thermohaline mixing prescriptions (Sec. 6 and Sec. 7). Our results and their implications are discussed in Sections 7 and 8.

The following section presents a derivation of the fluid parameter known as the "reduced density ratio," which is a fundamental parameter in fluid models of thermohaline mixing, and measures the stability of a given fluid configuration to the instability that drives thermohaline mixing. Section 3 goes into a similar level of detail regarding different prescriptions for thermohaline mixing in 1D stellar models to the end of summarizing the prescriptions currently implemented in MESA, and demonstrating key differences between these prescriptions. As was noted in Sec. 1.2, while the amount of chemical mixing predicted by these prescriptions does depend on the density ratio (or, equivalently, the reduced density ratio), we stress that the density ratio characterizes stability and thus the tendency to mix, not the actual amount of mixing or mixing efficiency (which is characterized by D th , introduced in Sec. 3).

The instability driving thermohaline mixing requires a Ledoux-stable inversion of the mean molecular weight µ stratification in the presence of a stable temperature gradient. The stability of the temperature gradient is given by the Schwarzschild criterion:

where the temperature gradient ∇ ≡ d ln P/d ln T (pressure P and temperature T ) has an adiabatic value ∇ = ∇ ad and saturates to ∇ = ∇ rad in hydrostatically stable regions where the flux is carried radiatively. The Ledoux criterion for convective stability is (Ledoux 1947)

which must be satisfied despite the inversion of the mean molecular weight. The Ledoux criterion accounts for the composition gradient ∇ µ = d ln µ/d ln P , where δ = -(∂ ln ρ/∂ ln T ) P,µ and φ = (∂ ln ρ/∂ ln µ) P,T (where ρ is density). To be unstable, the destabilizing µ gradient (the free energy source that drives thermohaline mixing) must be sufficiently strong relative to the stabilizing influence of the temperature gradient. This is quantified by the density ratio,

where R 0 < 1 implies the µ gradient is sufficiently unstable to drive convection, and R 0 > 1 implies the fluid is stably-stratified (i.e. no convection)foot_0 . As reviewed by Garaud (2018), fluids with R 0 > 1 can be prone to only if ∇ = ∇ rad . Thermohaline mixing primarily mixes chemicals, but does produce some minimal thermal mixing (see, e.g., Fig. 4 of Brown et al. 2013); thus, ∇ = ∇ rad . This thermal mixing is often ignored in mixing prescriptions in 1D stellar evolution programs, however.

double-diffusive instabilities whenever the thermal diffusivity, κ T , is greater than the compositional diffusivity, κ µ . Specifically, the instability driving thermohaline mixing acts whenever (Baines & Gill 1969) where

Note that in stellar radiation zones, typically τ 10 -6 . This means that very slight inversions of µ (large R 0 ) can drive thermohaline mixing, even when the temperature gradient is strongly stable according to the Schwarzschild criterion.

Throughout this paper, we express the density ratio R 0 in terms of the reduced density ratio

per, e.g., Traxler et al. (2011); Brown et al. (2013). In terms of r, the condition for thermohaline instability, Eq. ( 5), is

where r ≤ 0 is the threshold for convection and r ≥ 1 corresponds to scenarios where the µ inversion is too weak to drive the thermohaline instability.

Equation (8) can be readily evaluated at any radial location in a model star generated with a 1D stellar structure and evolution program. However, predicting the efficiency of thermohaline mixing is much more challenging. A diffusive approximation is commonly taken for the turbulent mixing of chemicals such that the total mixing of chemical species is given by the sum of the molecular diffusivity and a turbulent mixing coefficient, D th . This diffusion coefficient, D th , relates the rate of chemical mixing to the chemical composition gradient, and is broadly what is meant when referring to "mixing efficiency" throughout this paper. This quantity can be converted to the compositional Nusselt number discussed in the fluid dynamics literature, Nu µ , using the formula

We call any model that predicts D th as a function of stellar structure variables (e.g. gradients and molecular diffusivities of chemicals and heat) a parameterized mixing model or mixing prescription. Efforts to develop thermohaline mixing prescriptions for use in mod-els of stellar interiors date back many decades, see Garaud (2018) for a full review. Such mixing prescriptions have been implemented in a variety of 1D stellar evolution programs (see Lattanzio et al. 2015, and references therein), enabling studies of the effects of thermohaline mixing in stars across the Hertzprung-Russell diagram. Here, we briefly summarize the most commonly used and more recently developed prescriptions.

The de facto thermohaline mixing model used in MESA (first described in Cantiello & Langer 2010 and implemented for public use in Paxton et al. 2013) is commonly referred to as the "Kippenhahn model" and was originally derived by Ulrich (1972) and Kippenhahn et al. (1980). Using arguments based on dimensional grounds and assumptions about the shapes and motions of discrete fluid parcels, they derived a mixing efficiency of the form

(see Eq. ( 5) of Charbonnel & Zahn 2007) where C t is a free parameter, with plausible values ranging from C t = 658 (Ulrich 1972) to C t = 12 (Kippenhahn et al. 1980). We note that Eq. ( 10) predicts finite mixing for r ≥ 1 (R 0 ≥ 1/τ ), even though thermohaline mixing is formally stabilized for these parameters. Nevertheless, Eq. ( 10) is implemented in MESA as

(see Eq. ( 14) of Paxton et al. 2013). Here, α th is a dimensionless efficiency parameter related to C t by C t = 3α th /2, K is the radiative conductivity, ρ is the density, and C P is the specific heat at constant pressure, with κ T = K/(ρC P ). The green curve in Fig. 1 shows D th /κ µ vs. r calculated according to Eq. ( 11) for τ = 10 -6 , which is a representative value for the thermohaline-unstable region of RGB stars, and the same α th = 2 assumed in Cantiello & Langer (2010).

In addition to tension regarding the choice of model parameters (e.g. α th ) controlling overall mixing efficiency within a given 1D prescription (see e.g. Ulrich 1972;Kippenhahn et al. 1980;Charbonnel & Zahn 2007;Cantiello & Langer 2010;Traxler et al. 2011, for further discussion), there have also been questions about the appropriate trends in mixing efficiency as a function of fluid parameters (i.e. how D th should depend on quantities like r and Pr) and therefore the stellar structure variables on which they depend (Garaud 2018). Thus, recent work has sought to refine these mixing prescriptions by performing numerical experiments with multidimensional simulations to more accurately parameterize mixing efficiency (Denissenkov 2010;Traxler et al. 2011). Traxler et al. (2011) and Brown et al. (2013) performed 3D hydrodynamic simulations across a broad range of parameters. Not only did they find orders of magnitude less mixing than what is predicted by the Kippenhahn model with the model parameter required in Charbonnel & Zahn (2007) to find agreement with observations (C t = 1000), they also developed new mixing prescriptions that fit their simulations much more closely. In the case of Traxler et al. (2011), the authors derived a mixing law by fitting an analytic function of the form

to their simulation results, where

is the Prandtl number, with ν the kinematic viscosity, and a, b, and c are constants which they fit to data. While Traxler et al. (2011) clearly showed their simulations are inconsistent with the mixing efficiency D th implied by the Kippenhahn model with α th , C t ∼ 10 2 -10 3 , it is important to note that their simulations generally explored Pr, τ ∼ 10 -1 , whereas these fluid parameters are generally of the order 10 -6 in the radiative interiors of RGB stars. Thus, a fair question is whether mixing efficiency might increase to these larger values as Pr and τ approach 10 -6 . However, Traxler et al. (2011) varied these parameters by an order of magnitude in their simulations, and investigated trends in D th . They found that mixing should not increase in this fashion, as indicated by the dependence of D th on Pr and τ in Eq. ( 12), which makes an argument that these models can be made to fit the observational data difficult to justify. Brown et al. (2013) note that the model in Eq. ( 12) performs well at high R 0 , but underestimates mixing at low R 0 , particularly in the stellar regime of low Pr and τ . They develop a semi-analytical model for thermohaline mixing,

where λ is the growth rate of the fastest-growing linearly unstable mode, l is its horizontal wavenumber, and C ≈ 7 was fit to data from 3D hydrodynamic simulations. Both λ and l are functions of Pr, τ , and R 0 , and are obtained by finding the roots of a cubic and quadratic polynomial (their Eqs. 19 and 20). The orange curve in Fig. 1 shows D th /κ µ vs. r calculated according to Eq. ( 14) for Pr = τ = 10 -6 , representative values for the thermohaline-unstable regions of RGB stars. Note that D th /κ µ → 0 as r → 1 as expected, since the thermohaline instability becomes stable for r ≥ 1. We see that Eq. ( 11) with α th = 2 agrees with this prescription for some values of r, suggesting that significantly larger values of α th are not consistent with 3D hydrodynamic simulations. While the general dependence of D th /κ µ on r is significantly different between these two models, they do both feature monotonically decreasing values of D th /κ µ versus r. This prescription is implemented in MESA and has since been used in Bauer & Bildsten (2019) and other works. Harrington & Garaud (2019) extended the work of Brown et al. (2013) by performing 3D magnetohydrodynamic (MHD) simulations of thermohaline mixing with initially uniform, vertical magnetic fields of various strengths to approximate the effects of magnetic fields from external processes including, for instance, a global dipole field or a large-scale magnetic field left behind by a dynamo acting in the receding convective envelope. They found that magnetism strictly increases mixing efficiency, sometimes dramatically. They developed a mixing prescription that accounts for this effect by building on the model of Brown et al. (2013). The strength of the magnetic field is introduced into their model through their parameter H B , which is proportional to the square of the magnetic field strength and depends on other stellar structure variables (see Eq. 19 of Harrington & Garaud 2019). Their mixing prescription is of the form

where w f is obtained by solving a quartic polynomial that includes the magnetic field strength through H B , and K B 1.24 is directly related to the constant C in Eq. ( 14).

This mixing prescription agreed remarkably well with their 3D simulations, which were limited to r = 0.05 but ranged in magnetic field strength over several orders of magnitude. The prescription, which has not yet been implemented in MESA at the time of this writing, has two asymptotic limits, one where ŵ2 f ∝ B 2 0 when the magnetic field strength B 0 is large, and one which reduces to the model of Brown et al. (2013) when B 0 is small.

The purple curve in Fig. 1 shows D th /κ µ vs. r calculated according to Eq. ( 15) for the same parameter choices as the orange curve, and with H B = 10 -6 , appropriate for the thermohaline zone of a 1.1 M star at [Fe/H] = -0.2 and a magnetic field whose strength is O(100 G). Note that this magnetic field strength dra-matically increases mixing efficiency relative to the hydrodynamic values, particularly at larger values of r, whereas the model predicts the same mixing as the Brown model for r 10 -5 . For larger values of r, the dependence of D th /κ µ on r is profoundly different than either of the hydrodynamic models, with D th /κ µ increasing with r, even as the thermohaline instability approaches marginal stability as r → 1.

The dramatic enhancement in mixing efficiency predicted by this model for magnetic field strengths of even O(100 G) presents a promising resolution to the tension discussed above, namely that 1D stellar evolution models can only reproduce observations by assuming that mixing is far more efficient than what is seen in 3D hydrodynamic simulations. However, while their prescription may predict mixing efficiencies that are comparable in overall magnitude to that of the Kippenhahn model with C t ∼ 10foot_2 , the two prescriptions yield qualitatively different trends in D th vs r. Given the variance of the predictions of these models, we focus in this paper on showing how observations can be used to suggest the trends in mixing that models should hope to explain rather than on trying to calibrate model parameters controlling the overall mixing efficiency for a particular model, as has been done before, in order to provide a framework in which we can distinguish between mixing prescriptions.

We use MESA stable release version 21.12.21 to conduct 1D numerical simulations of stars incorporating the effects of thermohaline mixing for metallicities ranging from [Fe/H] = -1.4 to 0.4 (Z = 0.00068 to 0.038) and masses from 0.9 to 1.7 M at resolutions of 0.2 dex and 0.2 M , respectively. We adopt the solar abundance scale of Grevesse & Sauval (1998) and the corresponding opacities of Iglesias & Rogers (1996), with low-temperature opacities of Ferguson et al. (2005). Our models assume a helium abundance and heliumenrichment ratio of Y = 0.2485 and dY dZ = 1.3426, respectively, as in Tayar et al. (2022). We use an Eddington T-τ relation for the atmospheric surface boundary conditions. We adopt the mixing length theory (MLT) prescription of Cox (1980) with a fixed value of α MLT = 1.6 times the pressure scale height (H p ). We use the Ledoux criterion for convective stability and neglect the effects of convective overshoot (Ledoux 1947). We use the pp_extras.net nuclear reaction network, which contains 12 isotopes. More details are available in Appendix A, and the exact configuration of our physical and numerical parameter choices is available on the associated GitHub repositoryfoot_1 . MESA history files for each of the 200 simulations discussed in Sec. 5 are available on Zenodo (Fraser et al. 2022).

Simulations are evolved at 1.25× the default mesh (structural) resolution and 2× the default time resolution on the pre-main sequence and main sequence. Once the models ascend the red giant branch and reach a surface gravity log g ≤ 3, resolutions are increased to 2× the default spatial resolution and 10× the default temporal resolution, respectively. Optimal resolution values were determined according to the convergence tests detailed in Appendix B.

We study four grids of stellar evolution simulations with different thermohaline mixing prescriptions. One grid employs the Brown et al. (2013) prescription 3 . while the other three employ the Kippenhahn et al. (1980) prescription with coefficients α th ∈ [0.1, 2, 700]. The Kippenhahn α th = 0.1 model has inefficient mixing and represents a regime where thermohaline mixing is present but weak. We study α th = 2 because (1) it is consistent with the default implementation in MESA and discussion in the instrument paper (Paxton et al. 2013); (2) it is used in previous work (Cantiello & Langer 2010;Tayar & Joyce 2022); and (3) it is consistent with findings from 2D and 3D hydrodynamical simulations in stellar regimes (Denissenkov et al. 2010;Traxler et al. 2011;Brown et al. 2013). We also study a more traditional α th = 700, which is of the order of literature values calibrated using stellar observations (Lattanzio et al. 2015;Charbonnel & Zahn 2007). We note that the timestep size of our α th = 700 cases violate the CFL timestep constraint detailed by Lattanzio et al. (2015); MESA implicitly timesteps diffusive terms and so this violation does not affect the stability of our solutions, and we discuss possible effects on accuracy in appendix B.

These three choices for α th in the Kippenhahn prescription also correspond to three different relationships between the timescale over which thermohaline mixing acts, t th , and the evolutionary timescale of the star, t evol . We take the thermohaline timescale to be the diffusive timescale associated with thermohaline mixing, t th = d 2 /D th , where d is the radial depth of the thermohaline zone. In the case where α th = 0.1, t evol t th , meaning the timescale for homogenization of the mixing zone is large. In the case where α th = 700, t evol t th and the homogenization timescale is short. In the intermediate case (α th = 2), the timescales are comparable; this is likewise true for the assumptions made in the Brown model, though their prescription does not involve α th . Given the range of mixing timescales probed, these models should confer some insight as to how (or whether) the inferred fluid parameters, including the reduced density ratio r, depend on both the input mixing timescale and the stellar parameters.

We wish to extract the reduced density ratio r. As described in Sections 1.2 and 2, r characterizes the stability of the mean molecular weight gradient in the thermohaline region above the burning shell. It measures the region's tendency to mix, not its mixing efficiency. We define a selection criterion that averages over many mass been introduced in the MESA implementation. Brown's model is reproduced by assigning the quantity thermohaline_coeff = 1, as was done here shells and many evolutionary time steps to ensure that our measured values of r are representative of thermohaline mixing during the relevant evolutionary regime.

To measure r in our MESA simulations, we first restrict to the appropriate evolutionary phase. We exclude all models for which MESA does not detect thermohaline mixing within m max ≤ m i < 1.1m max , where m i is the mass coordinate of the ith mass shell and m max is the mass coordinate coinciding with the instantaneous peak of the nuclear energy generation. The thermohaline zone extends from a maximum mass coordinate m heavy to a minimum mass coordinate m light with stratification ∆m = m heavy -m light . We exclude the first 21 models in which the thermohaline zone spans at least 10 mass shells. We then compute the evolution of ∆m of the jth model, δm j = ∆m j -∆m j-1 for model j and the previous 20 models and compute δm = (1/20) 0 j=-20 δm j to determine whether the mass of the thermohaline region has evolved appreciably over the past several timesteps. We expect δm to be relatively large when the thermohaline zone is developing and small when it is in a relatively steady state. We then measure = | δm /max(∆m j )|. If < 5 × 10 -3 , we consider the model to have reached a steady state of thermohaline mixing (classified as "good" or "stable"), at which point we compute r.

To compute the reduced density ratio r = (R 0 -1)/(τ -1 -1), we take the volume average r = r i dV i / dV i over a subset of mass bins i of the thermohaline zone. We volume-average the reduced density ratio r over the mass range bounded by m heavy + 0.1∆m < m i ≤ m heavy + (0.1 + 1/3)∆m. In the volume average, we set the volume element dV i = 4πr 2 i ∆r i and perform integration using the composite trapezoidal rule as implemented in NumPy (van der Walt et al. 2011). We stop extracting r after we have collected measurements over 1000 models, which captures the behavior of the saturated thermohaline zone and its eventual merging with the convective envelope. For each stellar evolution simulation, we report the median of the volume-averaged r over all of the stable models in which measurements were taken. Results are discussed in terms of the logarithm of this quantity, log 10 r.

A movie demonstrating the evolution of a thermohaline front and the reduced density ratio selection algorithm is available in Appendix C.

The reduced density ratio r is a measure which relates the star's stable thermal stratification to the destabilizing inverted mean molecular weight stratification in the thermohaline zone. Regions with small r (r < 1) are destabilized by a mean molecular weight inversion and will exhibit thermohaline mixing. Regions with large r (r ≥ 1) have a mean molecular weight inversion that is too weak to overcome the stable thermal stratification, and thus cannot drive thermohaline instability. Note that r alone does not determine the efficiency of thermohaline mixing. The efficiency is determined by another function, D th (r), such as the prescriptions laid out in Sec. 3. Most prescriptions show D th (r) decreasing with increasing r, but some new prescriptions (e.g., Harrington & Garaud 2019) suggest that D th (r) may increase with increasing r.

Our goal is to measure how r varies with stellar mass and metallicity in RGB thermohaline zones. Put differently, we want to understand how much nuclear reactions in the hydrogen burning shell destabilize the mean molecular weight gradient, and how large that instability is compared to the thermally stable background. Getting a robust measure of the average value of r in a thermohaline zone is difficult, because the thermohaline instability mixes the fluid, changing the mean molecular weight profile, and ultimately changing r in a manner according to the prescription D th (r) used. In order to understand both the "natural" value of r that the stellar model inherits (i.e. in the absence of significant mixing) and how choice of mixing models affect r, we run four grids of models. We run two grids where the thermohaline mixing timescale and the evolutionary timescale are comparable (Brown, Kippenhahn with α th = 2). We run a third grid where thermohaline mixing is fast compared to evolution (Kippenhahn, α th = 700), from which we can understand how rapid mixing affects the value of r. Finally, we run a fourth grid where evolution is fast compared to mixing (Kippenhahn, α th = 0.1), which allows us to probe r when mixing does not appreciably modify the composition profile.

Results from these four physical configurations are compared in Fig. 2: the upper left panel shows results from the Brown model; the remaining three show results from the Kippenhahn prescription with α th varying as indicated. The reduced density ratio log 10 r is shown as a function of mass and metallicity and indicated on the color bar and grid labels.

In all cases, the most notable trend is that log 10 r decreases along the diagonal from high masses and metallicities (upper left) to low masses and metallicities (lower right). There is particularly high qualitative similarity between the Brown model and Kippenhahn model with α th = 2, which correspond to similar thermohaline mixing timescales. The case with the lowest mixing parameterization is the Kippenhahn α th = 0.1 case, and there the span of log 10 r values is smallest. We also note that, unlike in the other three cases, log 10 r does not scale precisely monotonically with either mass or [Fe/H] in the Kippenhahn α th = 700 case. This makes sense, because this is the case where mixing is most efficient, and so r is measuring the results of the mixing prescription rather than e.g., the rate at which destabilizing 3 He is burned.

While there is no clear relationship between the spread of log 10 r values observed when using the Kippenhahn prescriptions and the values of α th adopted in each, there is a clear relationship between the median values of log 10 r and α th : the reduced density ratios are larger when mixing is highly efficient (i.e. t th t evol ). This is consistent with D th (r) in Eqn. 11; mixing increases r, which in turn decreases the mixing efficiency D th in the Kippenhahn model, and so an equilibrium between the destabilizing source and the mixing is reached at a higher value of r (and a corresponding lower efficiency in that prescription).

Most importantly, the overall behavior of log 10 r as a function of mass and [Fe/H] is consistent regardless of the thermohaline parameterization adopted. The robust trend across 1D thermohaline mixing model assumptions suggests that r may be useful to compare to physical data sets.

Note that r is a measurable parameter with a robust trend in metallicity and mass, r(M, [Fe/H]). In the next section, we will measure changes in [C/N] as a function of the same parameter space ∆[C/N](M, [Fe/H]), which is a proxy for D th (r) integrated over time. Then we will compare ∆[C/N] and r. A monotonic trend between these quantities suggests a simple monotonic relationship D th (r), and could be used to reduce the parameter space from a two-dimensional space (mass, metallicity) to a one-dimensional space (reduced density ratio, r). Furthermore, the direction of the trend can be used to validate thermohaline mixing prescriptions. For example, the purple HG19 line and the green Kippenhahn line in Fig. 1 would produce different trends in this comparison of observational mixing vs. r because of the opposite manner in which mixing efficiency depends on r.

As discussed in Section 3, we are trying to distinguish between different models of thermohaline mixing based not on the amount of mixing they predict, but on their trends as a function of the reduced density ratio r. As discussed in Section 5, this requires stars of a wide range of masses and metallicities. It is therefore quite convenient that modern spectroscopic surveys have recently begun collecting measurements of mixing diagnostics for large samples of stars whose masses are also well con- -3.73 -3.69 -3.55 -3.49 -3.80 -3.69 -3.60 -3.51 -3.45 -3.75 -3.67 -3.56 -3.48 -3.39 -3.70 -3.60 -3.48 -3.42 -3.34 -3.59 -3.51 -3.43 -3.33 -3.25 -3.54 -3.43 -3.35 -3.27 -3.18 -3.43 -3.32 -3.23 -3.16 -3.08 -3.34 -3.23 -3.14 -3.07 -2.97 -3.25 -3.14 -3.03 -2.95 -2.85 -3.13 -3.00 -2.91 -2.77 -2.71 0.9 1.1 1.3 1.5 1.7 -3.84 -3.78 -3.74 -3.71 -3.87 -3.80 -3.74 -3.71 -3.68 -3.84 -3.76 -3.71 -3.67 -3.65 -3.80 -3.71 -3.66 -3.64 -3.61 -3.76 -3.66 -3.61 -3.59 -3.56 -3.71 -3.62 -3.56 -3.54 -3.52 -3.64 -3.55 -3.50 -3.47 -3.44 -3.58 -3.50 -3.43 -3.40 -3.37 -3.51 -3.43 -3.36 -3.33 -3.29 -3.42 -3.35 -3.29 -3.24 -3.20 -3.74 -3.61 -3.54 -3.48 -3.76 -3.70 -3.57 -3.50 -3.43 -3.76 -3.67 -3.57 -3.47 -3.39 -3.70 -3.61 -3.47 -3.41 -3.35 -3.58 -3.50 -3.43 -3.35 -3.27 -3.51 -3.44 -3.35 -3.28 -3.22 -3.44 -3.35 -3.26 -3.19 -3.12 -3.37 -3.26 -3.19 -3.11 -3.03 -3.27 -3.17 Kipp, th = 700 -3.18 -3.26 -3.18 -3.14 -3.09 -3.14 -3.16 -3.12 -3. strained. We choose for this work to use the carbonto-nitrogen ratios, [C/N], measured from the Apache Point Galactic Evolution Experiment (APOGEE, Majewski et al. 2017). APOGEE is a Sloan Digital Sky Survey III and IV (Blanton et al. 2017) project using the 2.5 meter Sloan Telescope (Gunn et al. 2006) and the APOGEE spectrograph (Wilson et al. 2019) to obtain medium resolution (R ∼ 22,500) spectra of large numbers of stars across the galaxy (Zasowski et al. 2017;Beaton et al. 2021;Santana et al. 2021). These spectra are homogeneously reduced and analyzed using the ASPCAP pipeline (Nidever et al. 2015;Zamora et al. 2015;García Pérez et al. 2016) and the resulting stellar parameters are then calibrated using asteroseismic, cluster, and field data (Holtzman et al. 2015(Holtzman et al. , 2018;;Jönsson et al. 2020). We choose to use the APOGEE data because this calibration work has already been done and an asteroseismic overlap sample is already available to provide stars with precise accurate masses, though similar work could likely be done with, for example, the lithium abundances measured by the GALAH survey (Buder et al. 2019) or the 12 C/ 13 C data estimated from the APOGEE data using the Brussels Automatic Code for Characterizing High accUracy Spectra (BACCHUS, Masseron et al. 2016) pipeline (C. Hayes, submitted). We also note that the evolution of [C/N] is in some ways simpler for these low-mass stars than other mixing diagnostics. Unlike lithium, its abundance at the surface does not change significantly during the main sequence due to the effects of rotational and other mixing processes (Iben 1967). Its initial ratio seems to be somewhat metallicity dependent (Shetrone et al. 2019), with higher values at lower metallicity. As stars reach the first dredge up, there is a strong, rapid, mass-and metallicity-dependent change in the surface [C/N] ratio (Salaris et al. 2015;Masseron & Gilmore 2015;Martig et al. 2016;Ness et al. 2016;Spoo et al. 2022). The [C/N] ratio at the surface then remains constant until stars reach the red giant branch bump, after which there seems to be another rapid drop in the [C/N] ratio, particularly in stars of low metallicity (e.g. Gratton et al. 2000;Shetrone et al. 2019); it is this drop that has been associated with thermohaline mixing. For stars of particularly low metallicities, there are some suggestions of Upper RGB extra mixing (Shetrone et al. 2019), but this is not well motivated theoretically and is distinct from the processes we are discussing here.

To estimate the amount of extra mixing in these stars near the bump-which thermohaline models suggest should correlate with the mixing coefficient D th described above-Shetrone et al. ( 2019) estimated the drop in [C/N] just above the red giant branch bump. Their work used α-element enhanced, and therefore old and low-mass (∼0.9 M ), first ascent red giant branch stars and binned them in bins of 0.2 dex in metallicity. The location of the red giant branch bump was identified empirically as an overdensity of stars at a particular surface gravity in each bin. They then identified the log g regime around the red giant branch bump and fit a hyperbolic tangent function to measure the location and size of the drops in the [C/N] ratio. For simplicity, we have reproduced their results in Table 1.

We add to their analysis a sample of higher metallicity stars with asteroseismic masses from the APOGEE-Kepler overlap sample (APOKASC, Pinsonneault et al. 2014Pinsonneault et al. , 2018)). We do this because, according to our analysis in Section 5, higher mass, higher metallicity stars probe larger values of the reduced density ratio, r. We first bin the stars in mass (0.2 M ) and metallicity (0.2 dex). For consistency with Pinsonneault et al. (2018) and Shetrone et al. (2019), we use the Data Release 14 Table 1. Observed extra mixing drops in bins of mass and metallicity, corrected for the 0.1456 dex of unmixing observed that we assume is due to systematic errors. We also include the reduced density ratios calculated for each of these bins using the variety of models discussed in Section 5.

rKip,700 0.9 -1.4 ... 0.73 0.00013 0.00012 0.00015 0.00066 0.9 -1.2 ... 0.67 0.00016 0.00013 0.00017 0.00073 0.9 -1.0 ... 0.48 0.00018 0.00014 0.00017 0.00078 0.9 -0.8 0.52 0.36 0.00020 0.00016 0.00020 0.00099 0.9 -0.6 0.29 0.27 0.00026 0.00018 0.00026 0.00140 0.9 -0.4 0.16 0.21 0.00029 0.00019 0.00031 0.00143 We note however that while the abundance scale seems to shift between releases, the rank ordering does not change very much (Spoo et al. 2022), which means that the conclusions of this analysis are not strongly affected by the choice of Data Release or seismic parameters. Unlike in the Shetrone et al. (2019) analysis (e.g. their Figure 2), there is not a sufficient number of stars near the bump in each bin to detect and measure the extra mixing directly in the asteroseismic sample. Instead, we define a 'pre-mixing' bin of stars between log g of 3.4 and 2.8 dex whose oscillations have identified them as first ascent red giants (Elsworth et al. 2019), as well as a 'post-mixing' bin of RGB stars with surface gravities between 2.3 and 1.0 dex. We then compute the average [C/N] of stars in each of the pre-mixing and post-mixing bins. If both bins had at least three stars, then the difference between the pre-mixing and post-mixing average [C/N] is plotted in Figure 3. Because of the calibration and choices in the analysis pipeline (see e.g. Holtzman et al. 2018;Jönsson et al. 2020;Smith et al. 2021), the scale of the abundances, particularly for carbon and ni-trogen, is somewhat uncertain. There sometimes exist small trends with surface gravity and temperature that are not fully removed in the calibration process. This is notable in our measurement results here; in the highest mass, highest metallicity bins, we formally measure 'unmixing' near the red giant branch bump, i.e. an increase rather than a decrease in the average [C/N] near the red giant branch bump, which is inconsistent both with theoretical expectations and with measurements from other sources. Following discussions with the APOGEE team (C. Hayes, private communication), we have decided to correct for these effects by correcting the bin with the most 'unmixing' to have 0 mixing, and subtracting that change from all of the other bins under the assumption that the systematic measurement errors are consistent for stars of similar temperatures and gravities. Such calibrations are common in the literature (Holtzman et al. 2015;Buder et al. 2021), and because we are most interested here in the trend in mixing amounts as a function of the stellar parameters, we do not expect this shift to alter the results of this analysis, but we emphasize that care should be taken by future users of this data.

Given the measured amounts of mixing described in Section 6 and the reduced density ratios r computed for stars of various masses and metallicities described in Section 5, it is now possible to compare the observations to the predictions of various thermohaline models from Section 2 and to assess whether the observed mixing is qualitatively consistent with any such theoretical prescription. We show in Figure 4 the corrected changes in [C/N] compared to the inferred reduced density ratios on axes analogous to those of Figure 1, where the y axis represented the rate of mixing. The four panels correspond to the four modeling configurations described in Section 4.

We first note that the observed trends are not strongly sensitive to the assumed 1D mixing model: mixing decreases with increasing reduced density ratio r regardless of parameterization. Besides this, there are similarities and differences between the data and the theoretical predictions, which are reproduced in Figure 5 for comparison. A key finding is that the observed mixing is strongly correlated with the fluid parameter r as predicted; this is true for stars with different masses and metallicities but similar reduced density ratios. This parameter therefore reduces a 2D (mass, metallicity) parameter space to one. We observe a decrease in the amount of mixing as the density ratio increases, which is consistent with standard 1D prescriptions of thermohaline mixing but inconsistent with the prescription from Harrington & Garaud (2019), which was informed by magnetohydrodynamic simulations. We also find that the range of average reduced density ratios probed by the observational data we have available here is much smaller than the range of density ratios simulated by and studied within the theoretical 3D fluid dynamics community (e.g. Brown et al. 2013), although we note that the full range of ratios do appear in each individual simulation (see Appendix C).

Thermohaline mixing has long been considered the most likely candidate for explaining the evolution of the surface chemistry of low-mass upper red giant branch stars. In this analysis, we have shown that:

1. prescriptions informed by three-dimensional simulations of the thermohaline instability suggest that mixing rates should strongly depend on the reduced density ratio, r, but the shape of the correlation varies from model to model; 2. one-dimensional stellar evolution models suggest that the average reduced density ratio, r, should vary as a function of mass and metallicity, with stars of lower masses and lower metallicities having smaller r;

3. one-dimensional stellar evolution models suggest that the average reduced density ratio, r, is not strongly dependent on the assumed parameterized thermohaline mixing model (e.g. Brown vs Kippenhahn), and regardless of the choice of α th in the case of Kippenhahn;

4. one-dimensional stellar evolution models suggest that the average reduced density ratio, r, occupies only a small range (10 -4 < r < 10 -3 ) of the parameter space formally unstable to thermohaline mixing, 0 < r < 1;

5. observations suggest that the amount of mixing is strongly correlated with the reduced density ratio, r, in a way that qualitatively agrees with predictions from three dimensional simulations;

6. observations indicate that the mixing rate and reduced density ratio, r, are inversely correlated, a finding that is consistent with currently available 1D prescriptions informed by 3D simulations, but not consistent with the prescription put forth by Harrington & Garaud (2019) with fixed magnetic parameter H B .

Most importantly, we find that our proposed framework for connecting observations of extra mixing in red giants to the fluid simulations of the thermohaline instability through the medium of one dimensional stellar evolution models is feasible and robust. Whereas previous work has focused on calibrating individual thermohaline mixing prescriptions against observations, this framework allows for comparison between different mixing models themselves. It motivates more rigorous exploration into whether red giant branch extra mixing should be associated with the thermohaline instability, and it will facilitate the translation of observational information into theoretical simulations. As observational data sets improve, abundance measurements will only become more capable of constraining stellar interior mixing, the relevant fluid parameters, and potentially even the magnetic fields in these regions.

However, the present study is largely a proof of concept; there is room for significant development in all aspects of this method. From an observational perspective, we have so far only considered the change in [C/N] as an observational mixing diagnostic, and only in a fairly restricted set of stars. Looking at higher or lower masses, using a variety of mixing diagnostics, or probing stars including globular clusters, binary stars, and Observations are color coded by the metallicity bin of each data point. In general, there is a clear correlation between these parameters, suggesting that the observed mixing may indeed be related to the unstable mean molecular weight gradient. Mixing and the reduced density ratio r are inversely correlated, which is consistent with hydrodynamic thermohaline prescriptions.

dwarf galaxies-where the base composition and mixing history may be different-could all be illuminating expansions of this work. There is also the potential to look at the timing of the extra mixing, as well as mixing rate as a function of time. Further, it may be possible to investigate whether mixing can be connected to the sorts of asteroseismic diagnostics that probe the interior structure of a star, including its density profile (Kjeldsen & Bedding 1995), chemical discontinuities (Verma et al. 2017), internal rotation (Gehan et al. 2018), and magnetic fields (Bugnet et al. 2021).

On the modeling side, we have explored a coarse but reasonable range of possible thermohaline conventions that could have impacted our conclusions, but this parameter exploration was certainly not exhaustive. It is worthwhile and necessary to test whether different 1D modeling assumptions impact the direction of this trend, particularly in the case of other mixing-related physical assumptions. Key variations in this regard include the choice of prescription for the Mixing Length Theory (MLT) formalism and value for the mixing length parameter, α MLT , the treatment of convective boundary mixing and convective overshoot, and the choice of atmospheric surface boundary conditions-all of which are well known to affect thermodynamic quantities in the regime we study here (Tayar et al. 2017;Joyce & Chaboyer 2018a,b;Viani et al. 2018). For similar thermodynamic reasons, it is also important to explore more extreme metallicity regimes, as the global metal content dictates the behavior of ∇ µ . Further, the need to introduce rotational mixing alongside thermohaline mixing in 1D stellar models to achieve the desired observational reproductions is likewise well established in the literature (c.f. Charbonnel & Lagarde 2010), making the consider- The observed extra mixing near the red giant branch bump as a function of the reduced density ratio inferred from one dimensional stellar evolution models. While the conversion from the change in a mixing diagnostic to the fluid mixing rate is not trivial, and therefore we do not attempt it here, we note that the observed mixing amounts are strongly negatively correlated with r, with stars probing on average a relatively narrow range of the regime formally unstable to thermohaline mixing.

ation of rotational effects an obvious candidate for future theoretical work. Finally, while previous ongoing work has endeavored to find a theoretical justification for why thermohaline mixing should be so efficient that accounts for the entire amount of extra mixing observed after the RGB bump (by adding e.g. magnetic fields or rotation), this work sets a new target for theoretical efforts: in order for thermohaline mixing to be the primary source of extra mixing, not only should there be a mixing prescription that agrees with 3D simulations and also reproduces the amount of extra mixing, but such a prescription should also reproduce trends in this extra mixing as a function of fundamental stellar parameters like mass and metallicity. For instance, while Harrington & Garaud (2019) demonstrate that magnetic fields of moderate strength can dramatically enhance the efficiency of thermohaline mixing-plausibly to levels that explain the full amount of extra mixing observed after the RGB bump-this work shows that their mixing prescription does not predict trends in mixing rate as a function of the density ratio that are consistent with observations.

In short, this paper has demonstrated the viability of comparing observed signatures of extra mixing on the red giant branch to the predictions of models informed directly by fluid simulations, and it has done so in the framework of parameters used in those simulations. We therefore anticipate and welcome future explorations that will produce better understanding of whether extra mixing should indeed be associated with the thermohaline instability and how observations of real stars-which probe a fluid regime well outside the regime we are currently able simulate-can provide constraints on the physics of that instability.

The first four authors of this manuscript contributed equally to this work.

A. Fraser contributed the majority of text and was responsible for Secs. 2 and 3. M. Joyce wrote and tested the MESA modeling templates, including inlists, evolutionary and structural output, and run_star_extras functionality, wrote the data analysis and parameter extraction scripts in collaboration with E.H. Anders (Python), and contributed text. We note that Lattanzio et al. (2015) report in their Eqn. 23 the Courant-Friedrichs-Lewy (CFL) criterion for diffusive processes as the timestep criterion necessary to resolve the thermohaline process in stellar models. This timestep criterion is required for stability when turbulent diffusivities are explicitly timestepped, but is not required for stability when diffusion is implicitly timestepped, as it is in MESA (and see Dutykh 2016, for a more detailed discussion of CFLs and timestepping diffusive processes). Numerical stability does not imply accuracy; Tayar & Joyce (2022) studied the behavior of thermohaline mixing with α th = 2 in similar stars in MESA by decreasing the timestep from well above the CFL criterion limit to well below it, and found good qualitative agreement between solutions which violated the CFL criterion and those that did not. This gives us some faith in our results regardless of whether they are above or below this limit. However, we note that this timestep limit is inversely proportional to the turbulent diffusivity δt CFL ∝ D -1 th . For the Kippenhahn model (Eqn. 10) where D th ∝ α th , the CFL timestep limit is inversely proportional to the thermohaline mixing efficiency coefficient α th used. It is therefore computationally impractical for us to evolve models that do not violate the CFL criterion for the α th = 700 case; it took on the order of one day for Tayar & Joyce (2022) to evolve a model that did not violate the CFL criterion with α th = 2, so it would take on the order of one year for us to evolve a model with α th = 700. This is prohibitively expensive and beyond the scope of this work. We also note that MESA is designed for convergence rather than speed, and so is roughly 10 times slower than the Monash code used by Lattanzio et al. (2015), which contributes to the unfeasibility of these calculations (Cinquegrana et al. 2022).

We include the α th = 700 case in this work because it is a typical value used within the field, but we urge caution when using and interpreting the results of these simulations as they cannot be as rigorously tested as e.g., the Brown et al. (2013) model. We have confirmed by eye that both the highest mass, highest metallicity and lowest mass, lowest metallicity simulations with α th = 700 evolve in a logical manner: the thermohaline front propagates monotonically outwards and connects with the envelope CZ.

In Fig. 7, we plot the stellar structure vs. mass coordinate in the simulation which employs the Brown model and has M = 1.1M and [Fe/H] = -0.2. We limit the x-limits of the plot to the mass coordinate energy generation peak of the hydrogen burning shell on the left m max , and to 1.1m max on the right. An animated version of this figure is available online in the published HTML version of this article and on YouTube at https://youtu.be/XLU8aS2q5-o. 2013) mixing prescription. We zoom in near the hydrogen burning shell, as indicated by the pink eps_nuc line which shows the nuclear energy generation rate. The thermohaline region can be identified as the place where there is a large amount of mixing (the light green log_D_mix is large), and we shade the bulk of this region in grey. As described in section 4, we only measure r over 1/3 of this region by mass, and the region over which we take this measurement is shaded in a darker grey. Within this region, we plot r in dark green. A zoom-in on log10r within the dark grey region is shown in the inset, and the measured value of r identified by the algorithm is plotted as a dark green horizontal line. Additional plotted lines include R0 and 1/τ as described in Sec. 3. Various quantities are quoted in text at the top of the image, including the star's age in years, the thermohaline mixing timescale in seconds, the radial extent of the thermohaline zone in terms of both length (dr) and mass coordinate (dm), and a flag indicating that this model is identified as "good" by our algorithm. An animated version of this figure is available online in the published HTML version of this article and on YouTube at https://youtu.be/XLU8aS2q5-o.