A phenomenological modification of thermohaline mixing in globular cluster red giants

TL;DR

This paper proposes a modified thermohaline mixing scheme for low-mass red giants that better matches observed carbon and lithium abundances in NGC 6397, highlighting the need for further observational and theoretical work.

Contribution

It introduces a temperature-dependent diffusive scheme for thermohaline mixing that improves agreement with observed stellar surface abundances.

Findings

Modified diffusive scheme fits observed abundances in NGC 6397

Faster mixing in hotter regions, slower in cooler regions

Two stream mixing approach was less successful

Abstract

Thermohaline mixing is a favoured mechanism for the so-called "extra mixing" on the red giant branch of low-mass stars. The mixing is triggered by the molecular weight inversion created above the hydrogen shell during first dredge-up when helium-3 burns via 3He(3He,2p)4He. The standard 1D diffusive mixing scheme cannot simultaneously match carbon and lithium abundances to NGC 6397 red giants. We investigate two modifications to the standard scheme: (1) an advective two stream mixing algorithm, and (2) modifications to the standard 1D thermohaline mixing formalism. We cannot simultaneously match carbon and lithium abundances using our two stream mixing approach. However we develop a modified diffusive scheme with an explicit temperature dependence that can simultaneously fit carbon and lithium abundances to NGC 6397 stars. Our modified diffusive scheme induces mixing that is faster than the standard theory predicts in the hotter part of the thermohaline region and mixing that is slower in the cooler part. Our results infer that the extra mixing mechanism needs further investigation and more observations are required, particularly for stars in different clusters spanning a range in metallicity.

Figures

Tables

Table 1. Table 1 : Models tested in Case 1: Independently changing v 𝑣 v and l 𝑙 l in the two stream advective mixing scheme.

Case	$f_v$	$f_l$	Name
$V$	0.10	1	$V0.10$
	0.33	1	$V0.33$
	0.50	1	$V0.50$
	2.00	1	$V2.00$
	3.00	1	$V3.00$
	10.0	1	$V10.0$
$L$	1	0.10	$L0.10$
	1	0.33	$L0.33$
	1	0.50	$L0.50$
	1	2.00	$L2.00$
	1	3.00	$L3.00$
	1	10.0	$L10.0$

Table 2. Table 2 : Models tested in Case 2: Changing v 𝑣 v and l 𝑙 l to maintain constant D t subscript 𝐷 𝑡 D_{t} in the two stream advective mixing scheme.

Case	$f_v$	$f_l$	Name
$DV$	0.10	10.0	$DV0.10$
	0.33	3.00	$DV0.33$
	0.50	2.00	$DV0.50$
	2.00	0.50	$DV2.00$
	3.00	0.33	$DV3.00$
	10.0	0.10	$DV10.0$

Table 3. Table 3 : Models tested in Case 3: Changing f u subscript 𝑓 𝑢 f_{u} and f d subscript 𝑓 𝑑 f_{d} in the two stream advective mixing scheme.

Case	$f_d$	$f_v$,$f_l$	Name
$FD$	0.01	1,1	$FD0.01$
	0.50	1,1	$FD0.50$
	0.99	1,1	$FD0.99$

Equations

$$H_v=f_{OS}{H}_P.$$

$$D_{OS}=D_{conv}{e}^{-2z/H_v},$$

$$D_t=C_{t}K(\varphi /\delta )\frac{-{\nabla }_{\mu }}{({\nabla }_{\mathrm{ad}}-\nabla )}.$$

$$C_t=\frac{8}{3}{\pi }^2{\alpha }^2,$$

$$D\equiv \frac{1}{3}vl,$$

$${\tau }_{\mathrm{thm}}=\frac{l}{v},$$

$${\tau }_{\mathrm{thm}}\equiv \frac{l^2}{3D_t}.$$

$$l=\alpha H_P,$$

$$v=\frac{3D}{l}.$$

$$F_u=F_d=\frac{4\pi r^2\rho v}{2}.$$

$$4\pi r_k^2{\rho }_k{v}_{u,k}{f}_u=4\pi r_k^2{\rho }_k{v}_{d,k}{f}_d.$$

$$v_{u,k}{f}_u=v_{d,k}{f}_d=\frac{v}{2}.$$

$${\alpha }_{2i-1}=\frac{1}{2f_d{m}_{2i-1}}4\pi r_k^2{v}_{d,k}{\rho }_k,$$

$$b_k=\frac{1}{2}\left[4\pi (r_{k+1}{)}^2{v}_{d,k+1}{\rho }_{k+1}-4\pi (r_k{)}^2{v}_{d,k}{\rho }_k\right].$$

$$\begin{array}{rl}{\beta }_{2i}& =\frac{v_k}{f_u{l}_k}+b_k\\ {\beta }_{2i-1}& =\frac{v_k}{f_d{l}_k}\end{array}\}\text{if bk is negative,}$$

$$\begin{array}{rl}{\beta }_{2i}& =\frac{v_k}{f_u{l}_k}\\ {\beta }_{2i-1}& =\frac{v_k}{f_d{l}_k}-b_k\end{array}\}\text{if bk is positive.}$$

$$v_{\mathrm{new}}$$

$$l_{\mathrm{new}}$$

$$D_{\mathrm{new}}=\{\begin{array}{ll}i\times D_t& \text{if }T>T_{\mathrm{crit}},\\ o\times D_t& \text{if }T<T_{\mathrm{crit}}.\end{array}$$

Full text

††pubyear: 2017††pagerange: A phenomenological modification of thermohaline mixing in globular cluster red giants–A phenomenological modification of thermohaline mixing in globular cluster red giants

1 Introduction

Mixing in stars is a complicated subject and modelling mixing events proves to be a difficult task. It is perhaps not surprising that there is observational evidence for mixing due to processes not predicted by theory. This is particularly the case for post-main sequence stars ascending the Red Giant Branch (RGB).

As a low-mass ($\lesssim 2.5\text{M}_{\odot}$) star ascends the RGB, the energy generated by the hydrogen burning shell (H shell) is transported to the surface easily via the convective envelope. This drives further expansion and deepens the mass location of the base of the convective envelope. The envelope eventually makes contact with the ashes of main sequence hydrogen burning and mixes these ashes to the surface. This process is called first dredge-up (FDU) and is the first major mixing event in a stars life (Iben, 1964).

FDU leaves behind a hydrogen abundance discontinuity that becomes very important for subsequent evolution (Iben, 1964). This discontinuity is believed to inhibit any extra mixing between the H shell and envelope (Mestel, 1953; Kippenhahn & Weigert, 1990; Chanamé et al., 2005). As the H shell converts hydrogen into helium, the inert helium core grows in mass. During core growth the H shell will move out (in mass) until it eventually reaches the hydrogen discontinuity left behind by FDU. The amount of available fuel and opacity of the H shell increases and the rate of hydrogen burning drops, which causes a drop in luminosity but a rise in effective temperature. The star appears to retrace its path in the Hertzsprung-Russell diagram (HRD) during this structural change. After the H shell passes through the hydrogen discontinuity, the luminosity will grow with core mass and the star will move up the RGB again almost along its original path. More stars are observed in this region of the HRD, producing a bump in the luminosity function (LFB, also known as the “RGB bump” or simply “bump”), which is a well-observed property of globular clusters (Kippenhahn & Weigert, 1990; Riello et al., 2003). After the LFB, the hydrogen discontinuity is erased and it is thought that extra mixing can then occur (Lattanzio & Wood, 2003; Charbonnel & Zahn, 2007; Karakas & Lattanzio, 2014). Low-mass stars have sufficiently long RGB lifetimes to allow the H shell to connect to the convective envelope via extra mixing mechanisms. This connection allows abundance changes occurring below the convection zone to be observed on the surface.

Globular clusters are an ideal laboratory for studies of extra mixing in stars. This is because globular clusters contain hundreds of thousands of low-mass stars, and the globular cluster population covers a huge range in metallicity. Globular clusters are also known to harbour multiple populations (Gratton et al., 2004, 2012), where these populations are correlated with spreads in the abundances of light elements (e.g., O, Na, Cottrell & Da Costa, 1981). Carbon and lithium abundances are useful indicators of mixing episodes during stellar evolution because they burn at different temperatures, meaning they are independently sensitive to temperature changes in the star and therefore reflect the temperature and mixing history of stars. The abundances of carbon and lithium in particular trace the evolutionary history of a star within a globular cluster (Briley et al., 1990; Gilroy & Brown, 1991; Gratton et al., 2000; Shetrone, 2003; Smith & Martell, 2003; Weiss et al., 2004; Martell et al., 2008; Lind et al., 2009). This is especially evident for NGC 6397, where the lithium abundance has been shown to vary with luminosity and a noticeable drop in the abundance is evident after FDU (e.g., Lind et al., 2009). Furthermore, there are determinations of the carbon abundances for NGC 6397 (Briley et al., 1990), which makes this cluster the ideal place for us to start our investigation.

A popular theory that offers a solution to extra mixing is thermohaline mixing. First introduced by Ulrich (1972) and then advanced further by Kippenhahn et al. (1980) in a general stellar context, an environment suitable for thermohaline mixing occurs naturally in low- to intermediate-mass stars after FDU (Eggleton et al., 2008). During core hydrogen burning, 3He is efficiently destroyed in the core but is produced outside the core where lower temperatures halt destruction. After FDU, 3He is homogenised in the outer layers and it is in this way that the surface value of 3He in low-mass stars increases. Once the outer layers are lost to the interstellar medium during later stages of evolution, the amount of 3He released is significant. This gives rise to the 3He abundance problem discussed by many authors since the mid 1990’s (e.g., Hata et al., 1995; Olive et al., 1995; Eggleton et al., 2006), where the amount of 3He observed in the interstellar medium matches what is predicted from Big Bang nucleosynthesis but is much lower than predicted from enrichment by low-mass stellar models. This is discussed in more detail by Eggleton et al. (2006).

Therefore if Big Bang nucleosynthesis is correct, then the theory of low-mass stellar evolution must be incomplete and some process (or processes) must destroy 3He in low-mass stars. Helium-3 can undergo fusion via 3He(3He,2p)4He once the H shell has advanced to regions where temperatures are high enough (e.g., Denissenkov & Herwig, 2004). This fusion process creates more particles than it consumes and produces a mean molecular weight ($\mu$) inversion (Eggleton et al., 2006). Charbonnel & Zahn (2007) linked this phenomenon to thermohaline mixing. A cooler layer of higher $\mu$ material sits atop a hotter, lower $\mu$ layer undergoing 3He fusion. “Fingers” of material from the hotter 3He fusion layer penetrate into the layer above, mixing material towards the surface.

There are other explanations of extra mixing in the literature. Most 1D stellar models do not include rotation, though it is known that stars rotate and this fact has been a source of discussion surrounding the (currently) unexplained mixing seen above the LFB. In three papers, Charbonnel & Lagarde (2010), Lagarde et al. (2011), and Lagarde et al. (2012) constructed a grid of models that included both rotation and thermohaline mixing, and combined these effects by adding the diffusion coefficients of each. This method does not account for any possible interaction between the two instabilities. This is discussed by Cantiello & Langer (2010) who state that the diffusion coefficient for rotationally-induced instabilities is much lower than the value for thermohaline mixing by several orders of magnitude. Canuto (1999) show that horizontal forces caused by rotation decrease the magnitude and efficiency of thermohaline mixing, though Medrano et al. (2014) found the opposite: the horizontal instabilities created by rotation increase the efficiency of thermohaline mixing in their 3D simulations. It is clear that there is much that is not understood in this area. Maeder et al. (2013) examine the interactions of instabilities (including thermohaline mixing) in rotating stars, an analysis that had been missing from the literature. Rotation alone cannot account for the extra mixing seen on the RGB (Chanamé et al., 2005; Palacios et al., 2006), yet a combination of instabilities (such as rotation, gravity waves, and thermohaline mixing) could prove to be a candidate (Denissenkov et al., 2009; Denissenkov, 2012). This is discussed further in §5.

Additionally, the stellar environment has consequences for thermohaline mixing, which was discovered in an oceanographic context that is incomparable (beyond a general sense) to stellar conditions. Cantiello & Langer (2010) show that the dimensionless Prandtl number, $\sigma$, the ratio of the kinematic viscosity to the thermal diffusivity, is very small in stars at 10*-6*. However, the Prandtl number in water (where thermohaline mixing was first discovered) is around 7 (Cantiello & Langer, 2010). Low Prandtl numbers (i.e. when $\sigma$ $\rightarrow$ 0) result in the (stellar) environment becoming unstable. This is shown in the simulations performed by Bascoul (2007) where fluid environments with low Prandtl numbers are dynamically unstable and turbulent. Whether thermohaline mixing can survive in such turbulent conditions is still under debate (Busse, 1978; Merryfield, 1995; Bascoul, 2007; Cantiello & Langer, 2010).

One-dimensional models of thermohaline mixing also find discrepancies with observations. Angelou et al. (2015) attempted to match their stellar models to observations of lithium and carbon in globular cluster RGB stars. They could not find a simultaneous match, which led them to the conclusion that the diffusive mixing scheme adopted for thermohaline mixing may not be complete or the most appropriate. However, important benefits of the thermohaline mechanism as a possible solution to extra mixing are that it naturally begins at the LFB and the thermohaline environment naturally occurs during normal RGB evolution. Indeed, it is also possible that the thermohaline instability triggers other instabilities that may then dominate mixing (this is discussed further in §5).

The aim of this study is to investigate models of transport by thermohaline mixing. In §2 we describe our numerical methods. In §3 we discuss the stellar models used. In §4 we state our test cases and show the results we obtain, and in §5 we discuss the results and state our conclusions.

2 Numerical Methods

2.1 The evolution code MONSTAR

All stellar models have been run using MONSTAR, the Monash version of the Mt. Stromlo stellar evolution code, as described in Angelou et al. (2015). High-temperature opacities are provided by the OPAL opacity tables (Iglesias & Rogers, 1996). Also included are fixed-metal distribution (OPAL type-1) tables that have the solar mixture composition of Grevesse & Noels (1993) and $\alpha$-element enhancement of Achim Weiss (Iglesias & Rogers, 1996). Tables variable in C and O content (OPAL type-2) are based upon abundances of Grevesse & Noels (1993). MONSTAR uses low-temperature (below $10^{4}$ K) tables with variable C and N content from Lederer & Aringer (2009) (Karakas et al., 2010; Angelou et al., 2015). The mass loss prescription on the RGB is from Reimers (1975) with $\eta_{R}=0.4$ (Constantino et al., 2016).

MONSTAR has a network of nine species (1H, 3He, 4He, 7Be, 7Li, 12C, 13C, 14N, 16O, and a tenth, inert, species). All key reaction rates included in MONSTAR, as well as more detail on the reaction rates used, can be found in Table 1 of Angelou et al. (2015). For the importance of timestepping and spatial resolution in MONSTAR (and evolution codes used by other groups) see Lattanzio et al. (2015).

Convection is modelled according to the mixing length theory (MLT; Böhm-Vitense, 1958) and we adopt $\alpha_{\rm MLT}=1.69$. Mixing of chemical species is modelled via a diffusion equation (Campbell & Lattanzio, 2008). The Schwarzschild criterion is adopted to define convective boundaries. Also included is overshoot beyond the formal boundary. In §2.1.1 we describe our convective model in more detail.

2.1.1 Convective overshoot

Convective overshoot is included following the procedure of Herwig et al. (1997), which inserts an exponential decay in the velocity of overshooting material. Let $H_{P}$ be the pressure scale height, $f_{OS}$ a scaling factor, and the velocity scale-height $H_{v}$ be

[TABLE]

We set $f_{OS}=0.14$ for material overshooting the convective boundary to match the observed location of FDU and the RGB bump for NGC 6397 RGB stars (discussed in more detail in §3.1.1). The diffusion coefficient for overshooting material, $D_{OS}$, at a distance $z$ from the convective boundary is

[TABLE]

where $D_{conv}$ is the diffusion coefficient at the convective boundary. Note that the only convective boundary in our models is at the base of the convective envelope.

The value of $f_{OS}$ in Equation 1 affects the depth of FDU and consequently the luminosity of the RGB bump. A high overshoot factor means that more material overshoots the convective envelope border and FDU is deeper. When the H shell advances in mass during normal RGB evolution, it will thus encounter the discontinuity caused by FDU sooner, resulting in a lower bump luminosity.

2.1.2 Thermohaline mixing in MONSTAR

The diffusion coefficient for thermohaline mixing is given by Ulrich (1972) and Kippenhahn et al. (1980) as

[TABLE]

The dimensionless parameter $C_{t}$ is formally given by

[TABLE]

where $\alpha$ is the aspect ratio of thermohaline “fingers”, $\phi=(\partial\ln\rho/\partial\ln\mu)_{P,T}$, $\delta=(\partial\ln\rho/\partial\ln T)_{P,\rho}$ ($\phi=\delta=1$ for an ideal gas), $K=4acT^{3}/(3\kappa\rho^{2}c_{P})$ is the thermal diffusivity, $a$ is the radiation density constant, $\kappa$ is the Rosseland mean opacity, $c_{P}$ is the specific heat capacity, and all other symbols have their usual meanings. The value of $\alpha$ is not known a priori, hence we treat $C_{t}$ as a free parameter.

2.1.3 Mixing timescales

In regions where mixing is modelled by a diffusion equation, a particle will travel with velocity $v$ (cm s*-1*) over a mixing length $l$ (cm). The diffusion coefficient is related to $v$ and $l$ by

[TABLE]

where $D$ is in cm2 s*-1*. We calculate the timescale $\tau_{\rm thm}$ of mixing in this region by using the relation

[TABLE]

or

[TABLE]

2.1.4 Nuclear burning

Lithium-7, which burns at around $2\times 10^{6}$ K, is a useful tracer of mixing and burning. This is especially the case for low-mass RGB stars that have extensive convective envelopes. The sudden drop in surface 7Li after FDU is one of the key observations that provides evidence for the existence of extra mixing (Charbonnel et al., 1998; Smiljanic et al., 2009; Lind et al., 2009). Beryllium-7, which captures an electron to form 7Li as part of the ppII chain, is an integral isotope in the Cameron-Fowler mechanism (a means for producing 7Li in stellar interiors; Cameron & Fowler, 1971). Beryllium-7 has a short half-life and is produced via alpha capture on 3He as the first reaction in the ppII and ppIII chains. Helium-3 itself is important because destruction of 3He via 3He(3He,2p)4He is the reaction that produces a local decrease in $\mu$, thus driving thermohaline mixing.

Another key piece of evidence for extra mixing is observations of [C/Fe]. This ratio also shows a sharp decline after the completion of FDU, confirming that there is a non-canonical mixing event occurring in the stellar interior (Suntzeff, 1981, 1989; Carbon et al., 1982; Trefzger et al., 1983; Langer et al., 1986; Briley et al., 1990; Gratton et al., 2000; Bellman et al., 2001; Shetrone, 2003; Martell et al., 2008; Shetrone et al., 2010).

2.2 The post-processing nucleosynthesis code MONSOON

Our stellar models have been constructed using MONSTAR and then processed by the Monash post-processing nucleosynthesis code MONSOON (Cannon, 1993; Campbell & Lattanzio, 2008; Karakas & Lattanzio, 2014; Karakas & Lugaro, 2016). MONSOON requires as input the radius, temperature, density, mixing length, and velocity as a function of interior mass and time. The reaction rates used in MONSOON are from the JINA REACLIB database (Cyburt et al., 2010).

Note that in each code we employ different mixing algorithms. In MONSTAR we employ a diffusive mixing scheme as outlined in §2.1.2. In MONSOON we employ an advective mixing scheme as outlined in §2.2.1.

2.2.1 Two stream advective mixing

Our two stream phenomenological advective mixing approach is shown in Fig. 1 (Cannon, 1993). We set the value of the diffusion coefficient $D_{t}$ in MONSTAR according to Equation 3 and calculate $v$ and $l$ for use in MONSOON using

[TABLE]

where $\alpha=\alpha_{\rm MLT}=1.69$ and is the same as used in MONSTAR, and from Equation 5

[TABLE]

The total mass flux $F$ across radius $r$ is zero in accordance with mass conservation. Using this information the mass fluxes in the up and down streams are calculated using

[TABLE]

Fig. 1 schematically shows the two streams in this mixing regime, with each stream being divided by horizontal levels $k$. The thermohaline mixing velocities according to Equation 9 are calculated at the boundary of each horizontal level. From the structure model, the density and radius are also known at the boundary of each level. To determine the velocity of material in each stream at each level $k$ ($v_{u,k}$ and $v_{d,k}$ for the up and down stream respectively), we denote the fractional cross-sectional areas of the downstream, $f_{d}$, and the upstream, $f_{u}$, for the thermohaline region (and all regions where mixing occurs). The values of $f_{d}$ and $f_{u}$ also satisfy $f_{d}+f_{u}=1$. To conserve mass it is required that up and down mass fluxes are equal, therefore the flux in each stream is

[TABLE]

The velocities $v_{u,k}$ and $v_{d,k}$ are determined from the thermohaline velocity in Equation 9 by

[TABLE]

The factor of $\nicefrac{{1}}{{2}}$ in the right-hand side of Equation 12 is a direct result of the mass flux equalling zero according to mass conservation. From this condition it is necessary that the flux in each stream must be equal as shown in Equations 10 and 11, therefore the total mass flux is split evenly between the two streams regardless of the cross-sectional area of each stream. A wider stream will consist of material travelling at a slower velocity and vice versa for a thin stream in accordance with Equation 12.

Showing the derivation for the downstream only (similar equations can be derived for the upstream), an equation for $\alpha$, the fraction of mass replaced in a downstream cell moving vertically per second on level $k$, is given by

[TABLE]

where $m_{2i-1}$ is the mass of cell $2i-1$.

The difference between the vertical mass flow into and out of a given cell must move horizontally if mass is to be conserved. This flow rate is denoted by $b$ and shown below for level $k$

[TABLE]

If $b=0$ we still allow for horizontal mass flow. This is calculated using standard mixing length concepts where a blob of material will travel a mixing length $l$ with velocity $v$ before losing its identity. Let $\beta$ be the fraction of mass replaced in each cell moving horizontally per second. The equation for $\beta$ must also take into consideration the mass flow according to $b$ (Equation 14). For each cell on level $k$, we define equations for the fraction of mass replaced in an upstream cell per second, say $\beta_{2i}$, and a downstream cell per second, say $\beta_{2i-1}$, in relation to $b$

[TABLE]

[TABLE]

Equations 15 and 16 state that if $b$ is negative, mass must flow from the upstream into the downstream. If $b$ is positive, mass must flow from the downstream into the upstream. We take $f_{u}=f_{d}=0.5$ except in §4.3.

3 Stellar models

Using MONSTAR and MONSOON we evolve our stellar models from the zero-age main sequence (ZAMS) to the helium flash, which terminates RGB evolution. To allow meaningful comparisons with observations we fit our stellar models to the globular cluster NGC 6397 because it is the only globular cluster to date that has observations for both carbon (Briley et al., 1990) and lithium (Lind et al., 2009) along the RGB. It is imperative to have observations of both elements because carbon and lithium burn at different temperatures and are therefore independent indicators of extra mixing, as stated in §1.

3.1 Fitting stellar parameters

We test two metallicities that span the range of observed metallicities for NGC 6397. We test [Fe/H] = $-2.03$ (Gratton et al., 2003) and $-1.82$ (Reid & Gizis, 1998). The metallicity given in Gratton et al. (2003) is very close to the Harris (1996, 2010 edition) value of $-2.02$. Lind et al. (2009) find [Fe/H] = $-2.10$, which is similar to results found by Castilho et al. (2000) ([Fe/H] = $-2.0\pm 0.05$), Gratton et al. (2001) ([Fe/H] = $-2.03\pm 0.02\pm 0.04$)111Where the first set of error bars are internal errors and the second set are systematic errors within the abundance scale of the 25 subdwarfs analysed by Gratton et al. (2001)., Korn et al. (2007) ([Fe/H] = $-2.12\pm 0.03$), and Husser et al. (2016) ([Fe/H] = $-2.120\pm 0.002$). The Reid & Gizis (1998) value of $-1.82$ is the upper bound of NGC 6397 published metallicities.

We use $[\alpha/\text{Fe}]=0.34$ according to Gratton et al. (2003). We determine $Z=0.000246$ using [Fe/H] = $-2.0$ according to Gratton et al. (2003) and $Z=0.00039$ using [Fe/H] = $-1.8$ according to Reid & Gizis (1998).

We construct stellar models with masses $0.79\text{M}_{\odot}$ and $0.8\text{M}_{\odot}$ and $Y=0.24$. These produce main sequence turn-off (MSTO) ages of 13.8 Gyr and 13.2 Gyr respectively, and RGB tip ages of 15.1 Gyr and 14.4 Gyr respectively with [Fe/H] = $-2.0$, which are good matches for the age of NGC 6397 found by Gratton et al. (2003) of $13.9\pm 1.1$ Gyr. A $0.8\text{M}_{\odot}$ stellar model with [Fe/H] = $-1.8$ produces a MSTO age of 13.3 Gyr and an RGB tip age of 14.6 Gyr.

3.1.1 FDU and RGB bump magnitudes

The FDU and RGB bump absolute visual magnitudes of NGC 6397 giants are $M_{V}=3.3$ and 0.16 respectively (Lind et al., 2009) using the bolometric corrections of Alonso et al. (1999). Lind et al. (2009) do not report uncertainties in their magnitude estimates. Nataf et al. (2013) find the V-band magnitude of the bump to be $12.533\pm 0.046$. We apply this uncertainty to the bump magnitude of Lind et al. (2009) and find the values of the FDU and bump magnitudes to be $3.3\pm 0.046$ and $0.16\pm 0.046$ respectively.

As stated in §2.1.1, the value of the overshoot factor $f_{OS}$ affects the predicted magnitude of the bump. We test a range of $f_{OS}$ values from 0.02 to 0.15 to match the magnitudes of FDU and RGB bump, with results shown in Fig. 2 for our $0.8\text{M}_{\odot}$, $Z=0.000246$, $\text{[Fe/H]}=-2.0$ stellar model.

We find the stellar models that best match the FDU and bump magnitudes observed for NGC 6397 using our two chosen metallicities are a $0.79\text{M}_{\odot}$ star with $Z=0.000246$ using [Fe/H] = $-2.0$ (Gratton et al., 2003), $Y=0.24$ and $f_{OS}=0.14$ (henceforth our $0.79\text{M}_{\odot}$ model), and a $0.80\text{M}_{\odot}$ star with $Z=0.00039$ using [Fe/H] = $-1.8$ (Gratton et al., 2003), $Y=0.24$ and $f_{OS}=0.10$ (henceforth our $0.80\text{M}_{\odot}$ model).

3.1.2 Hertzsprung-Russell diagram

We plot the giant branch of our best-fit models according to our overshoot results (Fig. 2) on a HRD in Fig. 3 compared to the observational HRD of NGC 6397 (Lind et al., 2009).

In Fig. 3 we see that the model that best matches the giant branch of NGC 6397 is our $0.79\text{M}_{\odot}$ model. This is also the best fit for FDU and the LFB as discussed in §3.1.1. The difference in metallicities between our $0.79\text{M}_{\odot}$ model and $0.80\text{M}_{\odot}$ model (along with initial stellar mass) likely contributes to the HRD discrepancies in Fig. 3, though this does not change our results or conclusions.

To allow analysis of burning timescales and elemental changes in our model, we choose an epoch just after surface abundance changes due to thermohaline mixing have occurred. The epoch chosen is shown in Fig. 4. All analysis of internal stellar properties is done at this epoch.

The locations of FDU and thermohaline mixing (“Thm”) on the HRD are when abundance changes on the surface are observed due to these events. Interestingly we see from the inset of Fig. 4 that abundance changes due to thermohaline mixing first appear on the surface at a lower luminosity than the LFB. Generally it is thought that thermohaline mixing coincides with the LFB on the HRD but this is not necessarily the case according to Fig. 4 and is also found by other groups222As shown in Fig. 4 the difference in $\log(\text{L}/\text{L}_{\odot})$ between thermohaline mixing and the LFB is only $\sim 0.01$. Observationally this is extremely difficult to distinguish, hence it is sufficient to approximate that the LFB and extra mixing occur at the same time. (Charbonnel & Lagarde, 2010). Extra mixing is governed by 3He fusion inverting the $\mu$ gradient whereas the LFB is due to the structural change in the star as the H shell passes through the hydrogen discontinuity caused by FDU. These events need not coincide with each other because they occur at different temperatures. The lithium abundance changes at lower temperatures before the H shell reaches the discontinuity. Hence abundance changes caused by thermohaline mixing occur slightly before the reversal in the HRD, which requires the H shell to have reached the discontinuity. There is also a delay between when thermohaline mixing begins and when it connects to the convective envelope as Fig 5 shows in more detail. The amount of time elapsed is dependent upon timestepping and spatial resolution in the stellar evolution code (for a detailed discussion of this, see Lattanzio et al., 2015). In Fig. 5 we can see that thermohaline mixing begins prior to the H shell connecting with the hydrogen abundance discontinuity.

3.2 Nuclear burning timescales in our standard model

The timescales discussed below in §3.2.1 have been determined using MONSTAR. We denote $D_{0}$ as the thermohaline coefficient for the standard case with $C_{t}=1000$ and test three variations of $D_{t}$ according to $D_{t}=f\times D_{0}$ where $f=\nicefrac{{1}}{{3}}$, 1, 3. We show the results in Figs. 6 and 7. The case of $f=3$ is mixing that is faster than the standard $f=1$ case, and $f=\nicefrac{{1}}{{3}}$ is mixing that is slower.

The mixing timescales for the three cases ($f$ = $\nicefrac{{1}}{{3}}$, 1, 3) and various nuclear burning timescales are shown in Fig. 6. Timescales are plotted as a function of radius. The nuclear burning timescales depend upon temperature and density, and remain the same regardless of $f$.

3.2.1 Lithium and carbon burning timescales in the thermohaline region

There are two limiting cases for the mixing: if it is infinitely fast then all abundances are homogeneous; if there is no mixing the shapes of the abundance profiles are due to burning alone (decreasing with increasing depth/temperature). Although not dominant at the temperatures in the thermohaline region333At a time just after the initiation of thermohaline mixing, temperatures range from $\sim 20$ MK at the base of the thermohaline region to $\sim 2$ MK at the base of the convective envelope (Fig. 8). as shown in Fig. 8, the ppII chain is still operating, destroying 3He by alpha capture via ${}^{3}\text{He}(\alpha,\gamma)^{7}\text{Be}$ and producing 7Be. As stated previously, faster mixing ($f=3$ in Figs. 6 and 7) produces a more homogeneous 3He profile than slower mixing ($f=\nicefrac{{1}}{{3}}$). The timescale for 3He alpha capture is longer than the mixing timescale regardless of $f$, producing a 3He abundance that does not vary dramatically with $f$ as shown in the middle panel of Fig. 7.

Alpha capture on 3He produces 7Be and this occurs predominantly in the interior where 3He burns fastest. A fast rate of mixing will bring 3He from the interior (where temperatures are higher) towards the surface faster. Therefore less 3He will burn, producing less 7Be. However, the 7Be will also be brought from the interior faster, and although it will be destroyed via electron capture at nearly the same rate at all temperatures, the amount of 7Be destroyed is less than when mixing is slower.

The rate of 7Be destruction is always faster than the rate of mixing for all tested $f$, therefore over the entire thermohaline zone the 7Be abundance does not become homogeneous. In the interior (from the base of the thermohaline zone to $\log(\text{R}/\text{R}_{\odot})\sim-0.8$), 7Li destruction is much faster than mixing by several orders of magnitude, therefore the 7Li profile is not homogeneous with position. The 7Li profile is complicated because both destruction (by proton capture) and production (by electron capture on 7Be) timescales are of similar orders of magnitude, particularly around $\log(\text{R}/\text{R}_{\odot})\sim-0.8$. This produces the variation of around 3 orders of magnitude that is dependent upon $f$ and shown in the top panel of Fig. 7. Beyond $\log(\text{R}/\text{R}_{\odot})\sim-0.8$ the rate of 7Li destruction is slower than the rate of mixing and the profile is approximately homogeneous.

Faster mixing means that the timescale of 12C destruction is slower than the mixing timescale over the entire thermohaline region and the abundance profile is (almost) homogeneous. When $f$ is reduced and the rate of mixing is slower, the timescale of 12C destruction can be faster than the mixing timescale for a thin region at the base of the thermohaline zone. More 12C destruction occurs because 12C is allowed to burn at the base of the thermohaline zone.

Destruction of 13C via proton capture is slower than both the rates of mixing and 12C destruction. When mixing is fast, 13C is essentially homogenised over the region. When mixing is slower, 13C is destroyed slower than it is produced, resulting in net production of 13C. The combination of increased 12C destruction and increased 13C production results in a 12C abundance that is lower and a 13C abundance that is higher in the thermohaline region when $f=\nicefrac{{1}}{{3}}$ compared to the $f=1$ and $f=3$ cases.

3.3 Motivation for studies and observational limitations

The data we use for NGC 6397 are from Lind et al. (2009) for lithium and Briley et al. (1990) for carbon. Angelou et al. (2015) found that the decrease in the carbon abundances due to thermohaline mixing occurred prior to the RGB bump whereas lithium declined at the bump, as predicted. We have found that this is related to the distance modulus and note that Lind et al. (2009) and Briley et al. (1990) used different values of the distance modulus for NGC 6397; Lind et al. (2009) report 12.57 and Briley et al. (1990) report 11.8. By adjusting the distance modulus of the stars observed by Briley et al. (1990) to match that of Lind et al. (2009), we resolve the issue found by Angelou et al. (2015) and find that carbon and lithium begin their decrease at the same magnitude.

To match lithium abundances to observations, a $C_{t}$ value of around 150 is required but a value closer to 1000 is needed to match carbon as shown in Fig. 9 (as also found by Angelou et al., 2015). To elaborate on this, if $C_{t}=1000$ is used and carbon matches observations then the models deplete lithium too much. Converse to this, if we use $C_{t}=150$ to match lithium observations then we require an increased depletion of carbon.

The spread of the observational data is a complication when comparing to our theoretical models (Fig. 9). The spread of the carbon abundances is of the order of 0.2 dex at a given magnitude. This spread is most likely due to the different stellar populations present in NGC 6397. Angelou et al. (2012) showed that observations of carbon and nitrogen in the intermediate-metallicity clusters M3 and M13 were best represented by separate models with initially different carbon and nitrogen compositions. The stars with an essentially normal composition are the CN-weak population, which form the upper envelope of the carbon distribution. It is these stars that we try to fit here. Errors of individual star carbon abundances are around 0.1 dex (Briley et al., 1990), as shown in the bottom panel of Fig. 9. The 0.1 dex error we use for the carbon observations is based upon the discussion and values given (in their Table 11) in Briley et al. (1990) and is a conservative value for the upper-RGB stars. In fact, the errors given in Briley et al. (1990) show their coolest stars towards the tip of the RGB have errors larger than their hotter stars, which makes sense given that the strength of the carbon molecular bands decreases with temperature. The model with $C_{t}=150$ can match the lower RGB stars within errors but cannot match the most carbon-depleted upper RGB stars as shown in Fig. 9.

The observations of lithium are more tightly constrained with smaller errors (Lind et al., 2009) but again there are very few observations near the RGB tip and the coolest few stars do not have associated errors444The observed abundances for the coolest RGB stars are upper limits only. as shown in the top panel of Fig. 9. When $C_{t}=150$ we can match the lithium abundances of all of the RGB stars in NGC 6397 but when $C_{t}=1000$ we cannot match any RGB stars after the start of extra mixing.

4 Test Cases and Results

4.1 Case 1: Independently changing $v$ and $l$ in the two stream advective mixing scheme

We explore the advective two stream model implemented in MONSOON by varying $v$ and $l$ in Equation 5. We invoke “mixing factors” to modify the mixing length and velocity according to

[TABLE]

where the subscript “std” indicates the parameter value for a given $C_{t}$ (as defined by Equation 4). Looking at the value of $\beta$ for horizontal mixing in Equations 15 and 16, we see that the only place in our two stream mixing model where $l$ appears is in the equations. Therefore it is straightforward to show that $\beta\propto 1/l$. Hence varying $f_{l}$ is the same as varying $1/\beta$.

It is important to note that when $f_{v}$ and $f_{l}$ are varied independently we change the effective value of $D_{t}$ in MONSOON by a factor of $f_{v}$ or $f_{l}$ despite setting a value of $D_{t}$ in MONSTAR according to the Ulrich/Kippenhahn formula. Therefore $D_{t}$ in MONSTAR and $D_{t}$ in MONSOON will be inconsistent. Indeed, as soon as we use a different mixing algorithm (in MONSOON) to the diffusion equation (in MONSTAR) then the results are technically inconsistent. We accept this inconsistency because there is negligible feedback on the underlying stellar structure.

Table 1 summarises the models calculated for Case 1. When referring to specific tests in Table 1 we first refer to the case letter and then the value of the variable. For example, when referencing the test where $f_{v}$ is being changed independently and we wish to refer to the $f_{v}=0.10$ test, we say $V0.10$.

4.1.1 Single parameter results

The effect of changing $f_{v}$ and $f_{l}$ independently on the surface lithium and carbon abundances is shown in Fig. 10 (for $f_{v}$) and Fig. 11 (for $f_{l}$).

Fast mixing and long mixing lengths (high $f_{v}$ and $f_{l}$ respectively) result in increased depletion of both carbon and lithium on the surface. Independently changing $f_{v}$ and $f_{l}$ cannot simultaneously match the lithium and the (upper envelope of the) carbon abundances to observations for a single value of $C_{t}$.

For the case of changing $f_{v}$ and $f_{l}$ independently when $C_{t}=1000$ the models that best match lithium abundances to observations are $V0.50$ or $L0.33$ (corresponding to “effective” $C_{t}$ values of 500 and 333 respectively). The models that best match carbon abundances to observations are $V1.00$ or $L1.00$ (corresponding to an “effective” $C_{t}$ value of 1000). These results are consistent with the motivation for this study (in §3.3) where if $C_{t}$ is chosen to match carbon to observations we deplete lithium too much even when uncertainties/errors of the observations are taken into account.

4.2 Case 2: Changing $v$ and $l$ to maintain constant $D_{t}$ in the two stream advective mixing scheme

Table 2 summarises the models tested for Case 2. Here both $f_{v}$ and $f_{l}$ are changed to maintain the selected value of $D_{t}$. We refer to the case letter and then the value of $f_{v}$ (remembering that $f_{l}=1/f_{v}$). For example, $DV0.10$ refers to the test where $f_{v}=0.10$ and $f_{l}=10$.

4.2.1 Results

Fig. 12 shows the effect on the surface abundances of lithium (top panel) and carbon (bottom panel). The effect of changing $f_{v}$ is much more significant than changing $f_{l}$. To elaborate, varying the velocity of the material has a greater effect on abundances than does varying the mixing length and this is why we see a similar result to changing $f_{v}$ independently (as in Fig. 10).

The model that best matches models to observations of lithium is $DV0.33$ and the model that best matches carbon is $DV1.00$. This is consistent with the results of independently changing $f_{v}$ and $f_{l}$ discussed in §4.1.1. In the single parameter tests it was found that models matched observations when $f_{v}$,$f_{l}\sim 0.33-0.5$ for lithium and when $f_{v}$,$f_{l}$ $\sim 1-2$ for carbon. We find the same result here. It is evident no solution can be found by modifying $f_{v}$ and $f_{l}$ to maintain $D_{t}$ using our two stream mixing algorithm.

4.3 Case 3: Changing $f_{u}$ and $f_{d}$ in the two stream advective mixing scheme

An important point to consider is that the cross-sectional areas of the streams are unknown; they can and may indeed be unequal in real stars (as opposed to $f_{u}=f_{d}=0.5$). The $f_{u}$ and $f_{d}$ values we tested are given in Table 3.

4.3.1 Results

Fig. 13 shows that modifying the up and down stream cross-sectional areas does not significantly affect the results. Indeed, only for the test $FD0.99$ do we see an effect on the lithium abundance, and only towards the end of RGB evolution.

We do no further tests of modifying $f_{v}$ and $f_{l}$ (either independently or to maintain $D_{t}$) for varying $f_{d}$ because Fig. 13 shows that modifying $f_{d}$ does not have a significant effect on the results and we expect the same result as seen in Figs. 10, 11, and 12 regardless of the value of $f_{d}$.

4.4 Case 4: Changing the thermohaline coefficient in the diffusive mixing scheme

The results from §4.1, §4.2 and §4.3 show that changing the mixing length/velocity in the advective scheme cannot match our models to observations. In the following we develop a method that can “target” (in a sense) different elements depending upon their destruction/production timescales.

Recalling that $D_{t}$ is related to $v$ and $l$ according to Equation 5, when we modify $D_{t}$ we effectively modify the mixing velocity and/or mixing length. Through the dimensionless value of $C_{t}$ we can produce different values of $D_{t}$ to match the observations, which could reveal characteristics of the extra mixing mechanism at the base of the convective envelope during RGB evolution. Manipulating $D_{t}$ as shown in §4.4.1 below expands the parameter space and increases the likelihood of being able to simultaneously match surface carbon and lithium abundances to observations. This may provide us with information on what is missing in the standard implementation. We perform these tests using MONSTAR with mixing modelled by a diffusion equation, as usual (described in §2.1.2).

4.4.1 Adding an additional temperature dependence

We can manipulate $D_{t}$ to modify abundances of 7Be, 7Li, 12C, and 13C in the thermohaline region by adding a new temperature dependence on $D_{t}$. We do this by setting a “critical temperature” $T_{\rm crit}$ with “inner” $i$ and “outer” $o$ factors such that

[TABLE]

From the base of the thermohaline region to the radial location of $T_{\rm crit}$ we multiply the diffusion coefficient $D_{t}$ as given in Equation 3 by a factor $i$. Similarly, from the location of $T_{\rm crit}$ to the envelope we multiply $D_{t}$ by a factor $o$.

The value of $T_{\rm crit}$ dictates the elements that are effected by $i$ and $o$. To elaborate, if a “high” $T_{\rm crit}$ value is selected so that it is sufficiently close to the base of the thermohaline region then $i$ will only affect elements burning at the highest temperatures and $o$ will have a larger effect on $D_{t}$ (and all abundances regardless of burning temperature). If $T_{\rm crit}$ is “low” and located sufficiently close to the envelope then $o$ will have little effect. If $T_{\rm crit}$ is somewhere in the middle of the thermohaline region then $i$ will predominantly affect the abundances of elements that have high burning temperatures (12C and 13C) and $o$ will predominantly affect the abundances of elements with low burning temperatures (7Be and 7Li).

4.4.2 Results

The number of combinations of $D_{t}$, $T_{\rm crit}$, $i$, and $o$ available mean that there are a family of solutions that can simultaneously match carbon and lithium abundances to observations. An analytic multi-dimensional theory of thermohaline mixing is not available and beyond the scope of this paper. Hence in the first instance we seek solutions that are within one order of magnitude of the diffusion coefficient as found by the Ulrich/Kippenhahn 1D theory. In Figs. 14 and 15 we show one solution where $T_{\rm crit}=8\text{MK}$, $i=3.0$ and $o=0.1$. We could certainly generate a better fit, but since the modification we used is purely phenomenological we feel that this would not add any insights. The present result tells us that to match both carbon and lithium simultaneously we need faster mixing in hot regions and slower mixing in cooler regions.

5 Discussion and Conclusions

We have modified the thermohaline mixing model to be able to match lithium and carbon abundances to observations of NGC 6397 red giants. We do this by adding an additional temperature dependence to the thermohaline diffusion coefficient. Our successful model proposes mixing that is faster than the standard 1D theory predicts near the base of the thermohaline region, and mixing that is slower further out in the region towards the base of the convective envelope.

The effect of our modification to the standard theory is to facilitate the Cameron-Fowler mechanism by reducing the decline of the lithium surface abundance. Our model achieves this by removing beryllium faster from the vicinity of the H shell and allowing it to remain longer in the outer part of the radiative zone (where it captures an electron to produce lithium). A similar effect can be obtained when one uses a constant diffusion coefficient, which was shown by Denissenkov & VandenBerg (2003).

Our successful method results in a number of non-unique solutions because of the size of the parameter space available. Despite this success, there are caveats that must be noted. Our model is 1D and subject to the uncertainties and limitations that are inherent in all stellar evolution modelling (e.g., Karakas & Lattanzio, 2014). Also, we are adding an additional temperature dependence that is not yet driven by physics. There may be physics in the 1D theory derived by Ulrich (1972) and Kippenhahn et al. (1980) that could drive such an extra temperature dependence that has not been identified. This is an avenue for further research and beyond the scope of this paper.

Modifying our advective scheme mixing length and velocity parameters independently (producing a change in $D_{t}$) and dependently (maintaining constant $D_{t}$) could not simultaneously match carbon and lithium for the same set of parameters. Our advective regime, unlike a diffusion equation, does not mix along a composition gradient and all elements present in the mixing regions are carried in the streams with the same velocity as the streams themselves. Fast mixing and long mixing lengths (high $f_{v}$ and $f_{l}$ respectively) result in increased depletion of both carbon and lithium on the surface because they, and all other elements in the mixing region, are brought down more quickly to high-temperature regions where they are burnt. Changing the velocity of the streams and/or the mixing length does not affect certain elements differently to others, therefore carbon and lithium are depleted more when the velocity and mixing length are increased (and vice versa when the mixing velocity and length are decreased). This is the limitation of this method and is the main reason why it is unsuccessful. Indeed, the standard diffusion implementation also fails because it shows these characteristics.

We look to results from multi-dimensional studies to inform 1D stellar modelling. The theory derived by Ulrich (1972) is one-dimensional. Further, it cannot adequately constrain the aspect ratio ($\alpha$ in Equation 4) of the fingers. Studies using 3D simulations of surface (Robinson et al., 2003; Steffen et al., 2005; Collet et al., 2007) and interior (Stancliffe et al., 2011; Ohlmann et al., 2017) convection zones in stars have found that the upstream velocity is slower and the cross-sectional area of the upstream flow is larger than the respective values for the downstream, i.e. $v_{u}<v_{d}$ and $f_{u}>f_{d}$. These effects do not appear in the standard theory.

Multi-dimensional simulations of thermohaline mixing have also been performed. The 2D simulations of Denissenkov (2010) show that fingers (corresponding to $\alpha\sim 7$) of material arise in the oceanic thermohaline environment but blobs with $\alpha\sim 0.5$ occur in the RGB case. Denissenkov (2010) achieved fingers of material (with $\alpha>1$) in their RGB case only for highly viscous environments (viscosities that are 4 orders of magnitude higher than in real RGB stars). Garaud & Brummell (2015) found shearing in their 2D simulations when the Prandtl number was less than 0.5 that was not seen in their 3D simulations, and concluded that for sufficiently low Prandtl numbers, 3D models are necessary to resolve the thermohaline environment.

Several other 3D simulations of thermohaline mixing have been conducted in recent years (Denissenkov & Merryfield, 2011; Traxler et al., 2011; Brown et al., 2013; Medrano et al., 2014; Garaud & Brummell, 2015). Denissenkov & Merryfield (2011) compared their 3D work to the 2D simulations of Denissenkov (2010) and confirmed the results from the 2D simulations, as well as finding the excitation of gravity waves in their oceanic case but not their RGB case. Secondary instabilities (e.g., gravity waves) triggered by thermohaline mixing have been found by other groups (Traxler et al., 2011; Garaud & Brummell, 2015). Three-dimensional simulations show that the shape of the thermohaline fingers changes as conditions come closer to representing real stellar conditions (i.e. as the Prandtl number decreases), with the thermohaline fingers becoming more like blobs (Traxler et al., 2011). Subsequent studies found that thermohaline fingers became blobs over time (Brown et al., 2013; Medrano et al., 2014; Garaud & Brummell, 2015).

In all of the multi-dimensional thermohaline mixing simulations above, the sizes of the down and upstreams were approximately equal to each other, and we found from our 1D study no significant effect when the stream size was altered. However, the simulated stellar environments are not entirely representative of the conditions in real stars. Real stars are likely to be much more turbulent because the Prandtl number is extremely small in reality ($\sim 10^{-6}$) compared to simulations, which typically have Prandtl numbers $\sim 0.1-0.01$ (Traxler et al., 2011; Brown et al., 2013; Garaud & Brummell, 2015). Additionally, the density ratio, the ratio of the (stabilising) entropy gradient to the (destabilising) compositional gradient, in simulations is generally $\sim 1.1$ and much lower than the RGB value of $1.7\times 10^{3}$ (Traxler et al., 2011; Brown et al., 2013; Garaud & Brummell, 2015). Brown et al. (2013) found that a density ratio of 1.1 resulted in larger, convective-like plumes, whereas a more turbulent environment with a density ratio of 3 was more finger-like, indicating that the density ratio in real RGB stars could produce fingers of material as opposed to blobs.

It is clear that we do not yet adequately understand the thermohaline mechanism. Additionally, to provide theoretical stellar models sufficient (observational) constraints, more observations of both carbon and lithium in globular cluster red giants are needed for stars in the same cluster (other than NGC 6397) covering a range in metallicity. This will help us to determine what the implementation we used to match models to NGC 6397 is telling us about necessary modifications to the standard theory.

Acknowledgements

K.H. acknowledges the financial support of the Australian Postgraduate Award scholarship and would like to thank the Stellar Interiors and Nucleosynthesis (SINs) group at Monash for their helpful discussions.

The authors thank the referee for their useful suggestions.