THERMOHALINE INSTABILITIES INSIDE STARS: A SYNTHETIC STUDY INCLUDING EXTERNAL TURBULENCE AND RADIATIVE LEVITATION

For many decades, the theory of thermohaline convection was developed in the framework of oceanography, as it represents one of the major physical factors leading to the large streams of the Earth's oceans (Stern 1960, 1967; Turner 1973; Gargett & Ruddick 2003). It was applied to astrophysics in only a very few papers, such as Ulrich (1972) and Kippenhahn et al. (1980, hereafter KRT80). At the present time, thermohaline convection is recognized as a major mixing process which occurs in stars in the presence of an unstable μ-gradient associated with a stable temperature gradient, provided that the μ-gradient is not large enough to drive dynamical convection. If a blob of stellar material begins to move down, it exchanges heat with its surroundings more rapidly than particles, and goes on falling, thereby triggering the instability (Vauclair 2004, and references therein).

Various situations lead to thermohaline convection in stellar interiors. They may be classified as follows.

• 1.  
  Accretion of heavy matter onto the star, which may come from planetary material (Vauclair 2004; Garaud 2011; Théado & Vauclair 2012) or from an evolved companion in a binary system (Stothers & Simon 1969; Stancliffe et al. 2007; Thompson et al. 2008).
• 2.  
  Local reduction of the mean molecular weight due to nuclear reactions, for example in red giants (Eggleton et al. 2006, 2008; Charbonnel & Zahn 2007; Charbonnel & Lagarde 2010; Stancliffe 2010; Denissenkov 2010).
• 3.  
  Heavy element accumulation induced by atomic diffusion (Théado et al. 2009) or helium accumulation induced by diffusion in a mass loss flux as in main-sequence helium-rich stars (Vauclair 1975).

The variations with depth of the radiative accelerations on specific elements can lead to their accumulation or depletion in various layers inside stars, which may have strong consequences on the stellar structure and evolution (Richer et al. 2000; Richard et al. 2001). The accumulations of iron and nickel, which represent important contributors to the opacity in some stellar layers, may lead to extra convective zones and can even, in some cases, trigger stellar pulsations through the iron-induced κ-mechanism (Charpinet et al. 1996; Pamyatnykh et al. 2004; Bourge & Alecian 2006). However, in all these previous studies, the thermohaline convection induced by heavy element accumulation was not taken into account. The first computations of this effect were given in Théado et al. (2009), who found that the accumulation of heavy elements was attenuated by the thermohaline mixing process, but not completely suppressed. The helium settling which occurs simultaneously with the heavy element accumulation induces a stabilizing contribution to the global μ-gradient, which helps to keep some accumulation.

Recent numerical simulations of thermohaline convection (Denissenkov 2010, hereafter D10; Traxler et al. 2011, hereafter TGS11) shed new light on these previous computations. As discussed by Théado & Vauclair (2012), the thermohaline mixing coefficients previously used in the literature have to be revised. Also, in the computations by Théado et al. (2009), the effect of the radiative acceleration on the falling blobs was not introduced. Finally, the influence on the thermohaline convection of an external turbulence, already addressed by Denissenkov & Pinsonneault (2008), is derived and introduced in the present context.

In the present paper, we analyze this situation and we give a new derivation of the thermohaline mixing coefficient, including radiative accelerations and/or external turbulence, consistent with the numerical simulations. This new coefficient will be used in the future to compute the evolution of the element abundances inside stars, including the simultaneous effects of atomic diffusion and thermohaline mixing. The introduction of thermohaline mixing in the computations may have important consequences on the stellar structure and evolution.

Thermohaline convection is not a real diffusion process. It is a special kind of convection which presents, at least in the case of oceans, elongated cells called "salt fingers." This convection, induced by an unstable chemical stratification, leads to chemical mixing. For this reason, and for simplicity, it is generally modeled in the same way as diffusion, with a mixing coefficient whose expression has to be deduced from physical considerations. In the past, the mixing coefficients proposed for stellar conditions were very uncertain, with unknown normalization values. For example, the expressions given by Ulrich (1972) and KRT80 differ by nearly two orders of magnitude. Recently, numerical simulations like those of TGS11 have led to real improvements, as the mixing coefficient must account for the chemical fluxes obtained from the simulations.

An important physical parameter in the study of thermohaline convection is the so-called density ratio in oceanography, which may be transposed to the astrophysical case as

$$\begin{equation} R_0 = \frac{\nabla _{{\rm ad}}-\nabla }{|\nabla _{\mu }|}, \end{equation} \tag{ 1 }$$

where ∇_ad and ∇ are the usual adiabatic and local (radiative) gradients dln T/dln P and ∇_μ = dln μ/dln P.

As discussed in Vauclair (2004), in the absence of external turbulence or radiative levitation, thermohaline convection may develop only if

$$\begin{equation} 1 < R_0 < \frac{1}{\tau }, \end{equation} \tag{ 2 }$$

where τ is the inverse Lewis number, equal to the ratio of the particle diffusivity to the thermal diffusivity.

These two boundaries, which are well known in oceanography, have important physical meanings. For thermohaline convection to occur, the medium must be stable according to the Ledoux criterion (no dynamical convection), but the relative variations of T and μ inside a falling blob must be large enough, compared to the surroundings, to trigger the instability. This last condition imposes that the ratio of the thermal diffusivity to the particle diffusivity, which induces the variations of these two quantities inside the blob, be larger than the external ratio of the temperature gradient to the particle gradient.

If R₀ is smaller than one, the μ-gradient triggers dynamical convection. If R₀ is larger than 1/τ, the initial μ-gradient is not large enough to overcome the stabilizing effect of the temperature gradient. This means that the expression for the thermohaline coefficient must behave in the right way at these two limits: it must become very large for R₀ close to one, and it must vanish for R₀ close to 1/τ. This behavior is not reproduced by any of the previously used thermohaline mixing coefficients, except for that proposed by D10, as discussed below.

Ulrich (1972) and KRT80 proposed to model the thermohaline mixing process in terms of a diffusion coefficient proportional to the inverse μ-gradient, which can simply be written (see Théado & Vauclair 2012) as

$$\begin{equation} D_{{\rm th}}=C_t k_T R_0^{-1}, \end{equation} \tag{ 3 }$$

where k_T is the thermal diffusivity and C_t is an unknown factor including the aspect ratio (length/width) of the fingers. C_t is evaluated as 12 by KRT80 and 658 by Ulrich (1972). In their computations of red giant stars, Charbonnel & Zahn (2007) considered this parameter as free and adjusted it to C_t = 1000. Note that this expression does not satisfy either of the two physical limits: D_th goes to infinity for R₀ = 0 and not R₀ = 1 as needed, and it does not vanish for R₀ = 1/τ.

The recent numerical simulations by TGS11 considerably improved our knowledge of thermohaline convection. They found universal scaling laws for the turbulent heat and compositional fluxes, valid for different sets of Prandt and Lewis numbers. Fitting their numerical results by an empirical law, TGS11 gave an expression for D_th which vanishes for R₀ = 1/τ, but reaches a finite limit for R₀ = 1, which is not correct. Interestingly enough, a close look at their simulation points (their Figure 2) shows that the heat and compositional fluxes do increase for small R₀, far above their empirical relation, suggesting that this relation should be revised.

D10 derived a thermohaline mixing coefficient using a straightforward linear treatment of Boussinesq equations in a quasi-compressible fluid. He found an expression which nicely satisfies the boundary conditions. Like Ulrich (1972) and KRT80, he introduced an adjustable parameter, C_t, depending on the shape of the convective cells (aspect ratio of the fingers) and tried to constrain it by two-dimensional simulations. Here we prefer to use the more recent three-dimensional simulations of TGS11 to adjust this parameter. We find that the D10 expression nicely fits the simulations, including the increase of D_th when R₀ becomes close to one, which is visible in the simulation but was not introduced in the TGS11 expression:

$$\begin{equation} D_{{\rm th}}=C_t k_T (R_0 - 1)^{-1}(1-R_0\tau). \end{equation} \tag{ 4 }$$

Figure 1 presents the behavior of all these coefficients as a function of R₀ in stellar conditions. As an example, they have been computed for a 1.7 M_☉, 403 Myr model. Here the 1/τ value, which is related to the rapid decrease of D_th for large values of R₀, is 6.8 × 10⁸, which is larger than the values generally given in the literature (e.g., TGS11 or Théado & Vauclair 2012). The basic reason for this difference is that this parameter is generally computed with helium as the diffusing element, whereas here we take into account the fact that the diffusive element is iron, which diffuses much more slowly than helium (cf. Figure 2). From all the considerations discussed above, the coefficient proposed by D10 is clearly the best one as it satisfies the results of the three-dimensional simulations as well as the boundary conditions. However, this coefficient does not take into account external turbulence, nor the effect of the radiative acceleration on the falling blobs, which may be important when the instability is induced by heavy element accumulation due to radiative levitation. In the following, we derive a new thermohaline coefficient in the same way as D10 but including these two effects.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this section, we give a linear analysis of thermohaline convection in the presence of atomic diffusion and/or external turbulence. We begin with Boussinesq equations that describe motions in a nearly incompressible stratified viscous fluid, in the same way as D10 except that we add into the equations the terms due to atomic diffusion and turbulent diffusivity.

Atomic diffusion is a selective process related to the fact that the stellar gas is composed of many different elements, each of them behaving in its own way in the stellar gravity and radiation fields. If the radiative acceleration g_R(i) acting on a given element (say iron) is locally larger than gravity, this element moves up with respect to the average stellar matter (composed mainly of protons). The element moves individually until it shares its acquired momentum with the surroundings, which occurs after a collision timescale. The overall effect of this individual radiative acceleration on the medium is of order g_R = Σ_iX_ig_R(i), where X_i is the mass fraction of the element considered. The global radiative acceleration on the medium g_R is generally several orders of magnitude smaller than the gravity, so that the medium as a whole is not destabilized (Figure 3). However, if the radiative acceleration on the element considered decreases with increasing radius, the diffusion flux decreases as well and this element accumulates in some layers inside the star. As a result, the internal stellar structure may be modified, due to changes in the opacity. Meanwhile, an inverse μ-gradient is created, which may become unstable against thermohaline convection.

Zoom In Zoom Out Reset image size

External turbulence (for example rotation-induced mixing) may also coexist with thermohaline convection. In previous works (e.g., Charbonnel & Lagarde 2010), such a situation was treated by simply adding the mixing coefficients of the two processes, computed separately, in spite of the analysis given by Denissenkov & Pinsonneault (2008), who showed that the thermohaline mixing coefficient must be modified in this case, due to the action of the external turbulence on the falling blobs. Here we recalculate and confirm this effect and introduce it into the final expression for the thermohaline mixing coefficient.

Following D10, we first write the Boussinesq equations in a nearly incompressible fluid in the following form:

$$\begin{equation} \frac{\partial \mathbf {v}}{\partial t}+(\mathbf {v},\nabla)\mathbf {v}=\frac{\delta \rho }{\rho _0}\mathbf {g_e}+\delta \mathbf {g_e}+\nu \nabla ^2 \mathbf {v} \end{equation} \tag{ 5 }$$

$$\begin{equation} \frac{\partial T}{\partial t}+(\mathbf {v},\nabla) T=k^*_T \nabla ^2 T \end{equation} \tag{ 6 }$$

$$\begin{equation} \frac{\partial \mu }{\partial t}+(\mathbf {v},\nabla) \mu =k^*_{\mu } [\nabla ^2 \mu + f(\mathbf {g_e}, T)]. \end{equation} \tag{ 7 }$$

In the first equation, we have to take into account the fact that the gravity acting on the fluid is now an effective gravity defined as g_e = g − g_R, where g_R is the overall radiative acceleration acting on the medium, namely g_R = Σ_iX_ig_R(i). As already mentioned, g_R is generally negligible compared to g. However, when the effective gravities inside and outside the blobs are subtracted, the g_R terms are the only ones remaining, as they are different in the falling matter and in the surroundings.

The second equation, which describes the heat transfer between the falling blobs and the external medium, remains unchanged except that, in case of external horizontal turbulence, the heat transfer induced by the extra mixing must be taken into account in the thermal diffusivity. When the Péclet number, the ratio of the horizontal turbulent diffusion coefficient to the thermal diffusivity, is much smaller than one, k*_T reduces to k_T, which is given by

$$\begin{equation} k_T=\frac{4 a c T^3}{3 \kappa c_P \rho ^2}. \end{equation} \tag{ 8 }$$

But when the Péclet number becomes larger than one, the horizontal diffusion coefficient must be added to k_T in the expression of k*_T.

The third equation concerns the diffusion of particles between the blobs and their surroundings. The coefficient k_μ is the so-called particle diffusivity which, in hydrodynamical contexts, does not take into account the fact that stars are made of multi-component gases. It is generally computed in an approximate way, and the values cited in hydrodynamical papers often correspond to the diffusion coefficient of helium in a hydrogen-rich medium. Sometimes it is even identified with the molecular viscosity, which is nothing else than the hydrogen–hydrogen self-diffusion coefficient. In the atomic diffusion context, the diffusion coefficients for each element are computed in a consistent way (e.g., Paquette et al. 1986). It is not possible, in the present situation, to take into account the separate diffusion of all the elements through the blobs' boundaries, but one must at least be careful to use for k_μ the diffusion coefficients corresponding to the most important elements (e.g., iron). The ratio of the diffusion coefficients for two different ions (1) and (2) of respective charges Z₁ and Z₂ is approximately (Z₂/Z₁)². Using that of iron leads to a smaller diffusivity than using that of helium, which would not be correct in this case (Figure 2).

In a general way, the particle diffusion between the interior and the exterior of the falling blobs is not only due to the abundance gradients, but also to the effects of gravity, radiative acceleration, temperature gradient, etc., which occur in the vertical direction (e.g., Vauclair & Vauclair 1982). This is added as $f(\mathbf {g_e}, T)$ in Equation (7). However, as pointed out by previous authors (e.g., D10), the vertical selective diffusion velocity is negligible compared to the downward velocity of the blobs so that we can ignore this effect. On the other hand, as mentioned above, the average effect of the radiative acceleration on the medium has to be taken into account in the computation of the velocity (Equation (5)).

Finally, in case of external turbulence, the mixing of particles between the blobs and the interblobs is increased. The horizontal component of the extra turbulent diffusivity has then to be added to k_μ in Equation (7). We have written the total diffusivity as k*_μ, which reduces to k_μ when no external turbulence is involved.

In the following, we use the Cartesian coordinate system with the vertical z axis pointing in the direction opposite to the gravitational acceleration, the horizontal x axis in a star's meridional plane and the horizontal y axis perpendicular to them. Linearizing Equations (5)–(7), with w representing the vertical component of the blob velocity, we obtain the following system of equations, similar to Equations (5)–(7) of D10 except for the gravity term and the effect of external turbulence included in k*_T and k*_μ:

$$\begin{equation} \frac{\partial w}{\partial t}=g_e(\alpha \delta T -\beta \delta \mu)+ \nu \nabla ^2 w + \delta g_e \end{equation} \tag{ 9 }$$

$$\begin{equation} \frac{\partial \delta T}{\partial t}=-w \frac{\partial T}{\partial z}+k^*_T \nabla ^2 \delta T \end{equation} \tag{ 10 }$$

$$\begin{equation} \frac{\partial \delta \mu }{\partial t}=-w \frac{\partial \mu }{\partial z}+k^*_{\mu } \nabla ^2 \delta \mu, \end{equation} \tag{ 11 }$$

where α and β are the usual coefficients of the equation of state:

$$\begin{equation} \frac{\partial \rho }{\rho _0}=- \alpha \partial T + \beta \partial \mu. \end{equation} \tag{ 12 }$$

Note that for simplicity we keep here the notation ∂T/∂z but we have to remember that in stars, due to the temperature stratification, this is equivalent to (∇_rad − ∇_ad)(T/H_P).

When a heavy blob falls down, the conditions inside it vary in a non-adiabatic way due to heat and particle diffusion with the surroundings. The vertical variations of T and μ inside and outside the blob may be expressed in the following way:

$$\begin{equation} \frac{\partial T}{\partial z}=\frac{\partial T_0}{\partial z}(1- \delta _T) \end{equation} \tag{ 13 }$$

$$\begin{equation} \frac{\partial \mu }{\partial z}=\frac{\partial \mu _{0}}{\partial z}(1- \delta _{\mu }), \end{equation} \tag{ 14 }$$

where the coefficients δ_T and δ_μ have values between 0 (equality between the internal and external values of T and μ) and 1 (adiabaticity, no internal variations of T and μ). They may be defined from the following relations which give the differences δT and δμ between the inside and the outside of the blobs:

$$\begin{equation} \delta T=- \frac{1}{2} \frac{\partial T_0}{\partial z} \ell \delta _T \end{equation} \tag{ 15 }$$

$$\begin{equation} \delta \mu =- \frac{1}{2} \frac{\partial \mu _{0}}{\partial z} \ell \delta _{\mu }, \end{equation} \tag{ 16 }$$

where T₀ and μ₀ are the temperature and mean molecular weight outside the blobs and ℓ is the unknown length of the elongated cells (fingers).

The blobs stop falling when the buoyancy inside and outside are identical, that is when ρg_e is the same inside and outside.

As in D10, we transform Equations (9), (10), and (11) in terms of w, δ_T, and δ_μ, and we search for solutions in the form of exp(σt) × exp(i[k_xx + k_yy]). But before that we have to express the new term δg_e in the same way (in terms of δ_T and δ_μ). In the expression for the effective gravity, the only term which is different inside and outside the blobs is the radiative acceleration, which is a function of both the temperature and the chemical composition. Although these two variables are not completely disconnected, we will separate them in a first approximation. This approximation will be discussed in the following section. We write

$$\begin{equation} \delta g_e= \frac{\partial g_e}{\partial T} \delta T + \frac{\partial g_e}{\partial \mu } \delta {\mu } \end{equation} \tag{ 17 }$$

or

$$\begin{eqnarray} \delta g_e &=& -\biggl (\frac{\partial g_e}{\partial T}\frac{\partial T_0}{\partial z} \frac{\ell }{2} \delta _T + \frac{\partial g_e}{\partial \mu } \frac{\partial \mu _0}{\partial z} \frac{\ell }{2} \delta _{\mu }\biggr)\nonumber\\ &=& -\biggl (\frac{\ell }{2} \frac{\partial \mu _{0}}{\partial z} \left[\frac{\partial g_e}{\partial T} R_0 \frac{T}{\mu }\delta _T + \frac{\partial g_e}{\partial \mu } \delta _{\mu }\right]\biggr). \end{eqnarray} \tag{ 18 }$$

In the following, for simplicity we assume a perfect gas law, so that α = 1/T and β = 1/μ. Also, for simplicity we write here k_T and k_μ instead of k*_T and k*_μ. Defining g_{e, T} as ∂ln g_e/∂ln T and g_{e, μ} as ∂ln g_e/∂ln μ, the preceding equations may then be reduced to

$$\begin{eqnarray} 4 \sigma ^3 &+& 2(k_T+k_{\mu }+2 \nu)k^2 \sigma ^2 \nonumber \\ &+& \left\lbrace [k_T k_{\mu }+2 \nu (k_T+k_{\mu })]k^4 + 2 g_e \frac{\partial ln\mu _0}{\partial z}(R_0-1)\right\rbrace \sigma \nonumber \\ &-& g_e \frac{\partial ln\mu _0}{\partial z}(k_T-R_0 k_{\mu }) k^2 +\nu k_T k_{\mu } k^6 \nonumber \\ &+& 2\sigma g_e \frac{\partial ln\mu _0}{\partial z}(g_{e,\mu } + R_0 g_{e,T}) + k^2 g_e \frac{\partial ln\mu _0}{\partial z}\nonumber\\ && \times (k_T g_{e,\mu } + R_0 k_{\mu } g_{e,T}) =0 \end{eqnarray} \tag{ 19 }$$

with k² = k²_x + k²_y.

Taking into account that the thermal diffusivity is much larger than the particle diffusivity and the molecular viscosity, we neglect the same terms as D10 and this equation is simplified to

$$\begin{equation} \sigma = \frac{k^2}{2} \frac{(k_T - R_0 k_{\mu }) + (k_T g_{e,\mu } + R_0 k_{\mu }g_{e,T})}{(R_0 - 1) - (g_{e,\mu } + R_0 g_{e,T})}. \end{equation} \tag{ 20 }$$

The thermohaline mixing coefficient, of the order of 1/2σℓ², becomes

$$\begin{equation} D_{{\rm th}}=C_t k^*_T \frac{(1 - R_0 \tau ^*) + (g_{e,\mu } + R_0 \tau ^* g_{e,T})}{(R_0 - 1) - (g_{e,\mu } + R_0 g_{e,T})}, \end{equation} \tag{ 21 }$$

where k_T has been replaced by k*_T and the inverse Lewis number τ has been replaced by τ* = k*_μ/k*_T. C_t is the usual constant related to the aspect ratio of the fingers, that we adjust using TGS11 simulations.

If the derivatives of the effective gravity are neglected, as well as the horizontal turbulence, this expression for the thermohaline mixing coefficient is identical to that of D10. If only horizontal turbulence is introduced, this expression is similar to that obtained by Denissenkov & Pinsonneault (2008).

In the general case, the two limits for the development of thermohaline convection are slightly modified. The boundary with dynamical convection, which occurs for large μ-gradients, is not exactly R₀ = 1 but R₀ = (1 + g_{e, μ})/(1 − g_{e, T}). The other boundary, for small μ-gradients, becomes R₀ = (1 + g_{e, μ})/[τ*(1 − g_{e, T})] instead of R₀ = 1/τ. The orders of magnitude of these corrections are discussed in the next section.

In the previous section, we have introduced the partial derivatives of the effective gravity g_{e, μ} = ∂ln g_e/∂ln μ and g_{e, T} = ∂ln g_e/∂ln T. We now have to understand their real meaning, and derive their orders of magnitude in stellar situations. In the following we first assume no external turbulence. We will discuss the effect of adding turbulence later on.

The effective gravity is g_e = g − g_R where g_R is the radiative acceleration transferred to the whole medium by the elements which are individually levitated, namely g_R = Σ_iX_ig_R(i) where X_i is the mass fraction of the element i and g_R(i) the radiative acceleration acting on the same element i. In the following, to evaluate the orders of magnitude of these terms, we will simplify the problem by assuming that one major element is the reason for the μ-gradient inversion (e.g., iron). Generalization to several elements will be done numerically in stellar evolution codes.

Even in the case of one levitating element, the situation is not simple as each ion behaves in a different way, which is precisely the reason for the selective accumulation. In this case, we write g_R = X_ig_R(i) where g_R(i) represents the contribution to the radiative acceleration of all the ions of the element considered, that is:

$$\begin{equation} g_R(i) = \Sigma _j X_{i,j} g_{i,j} = X_i \Sigma _j \frac{X_{i,j}}{X_i} g_{i,j}, \end{equation} \tag{ 22 }$$

where j represents the various ionization stages of the element i.

The variation of g_R with μ is a combination of two effects. First, the relative mass fraction of the heavy element X_i appears as a multiplicative coefficient in g_R. Second, this mass fraction has also an influence on g_R(i) through the saturation of the lines. Namely, when the abundance of the element increases, the effect of the saturated lines on the total radiative acceleration of the element considered decreases because all the ions considered have to share the same available photon flux.

The variation of μ induced by the accumulation of element i may be written:

$$\begin{equation} \delta \ln \mu \approx \delta {X_i} \left[1 - \frac{\mu \bar{n_i}}{A_i}\right], \end{equation} \tag{ 23 }$$

where $\bar{n_i}$ is the average number of particles (including electrons) associated with element i in the prevailing physical conditions and A_i is its mass number.

For the heavy elements considered, the second term is small compared to one (otherwise there would not be any inverse μ-gradient). We neglect it in a first approximation and write g_{e, μ} as

$$\begin{equation} g_{e,\mu } = \frac{\partial \ln g_e}{\partial \ln \mu } \approx \frac{1}{g} \frac{\partial g_R}{\partial X_i} = \frac{g_R(i)}{g} + \frac{X_i}{g} \frac{\partial g_R(i)}{\partial X_i}. \end{equation} \tag{ 24 }$$

In this equation, the first term is positive. The second term is negative and difficult to evaluate, but its maximum importance is reached if g_R(i) varies like 1/X_i, in which case g_{e, μ} vanishes. On the contrary, if the saturation effects are neglected, the g_{e, μ} term can reach the order of magnitude of g_R(i)/g, which may become larger than one (Figure 3).

The variation of g_R with temperature only appears through g_R(i) and comes mostly from the fact that the various ions of the same elements are not pushed up in the same way by the radiative acceleration due to their different atomic structure:

$$\begin{equation} g_{e,T} = \frac{\partial \ln g_e}{\partial \ln T} = \frac{T}{g} \frac{\partial g_R}{\partial T} = \frac{T}{g} X_i \frac{\partial g_R(i)}{\partial T} \end{equation} \tag{ 25 }$$

or:

$$\begin{equation} g_{e,T} = \frac{T}{g} X_i \Sigma _j \biggl (X_{i,j} \frac{\partial g_{i,j}}{\partial T} + g_{i,j}\frac{\partial X_{i,j}}{\partial T}\biggr). \end{equation} \tag{ 26 }$$

The order of magnitude of the first term may be evaluated with the assumption that g_{i, j} varies inside the star approximately like 1/T (Michaud et al. 1976). It is negative, of order −g_R/g. The second term is more difficult to evaluate, as X_{i, j} depends on the Saha equation. It may either increase or decrease with T. In the approximation of two ionization stages, j and (j − 1), the fraction of ion j is approximately proportional to exp (− E_{j − 1}/k_BT), where E_{j − 1} is the ionization energy of ion (j − 1), so that its contribution to the partial derivative of g_e may be approximately evaluated as (g_R/g)(E_{j − 1}/k_BT), which may become larger than g_R/g. In our reference model (1.7 M_☉ at 403 Myr), this term never exceeds 15% whereas g_R/g always remains below 3%.

It is interesting to note that, due to the signs of these terms, the effect of the radiative acceleration on the thermohaline mixing coefficient is often to reduce it, but may also sometimes increase it. In the regions where the radiative acceleration decreases with increasing radius, the falling blobs suffer less levitation than the surroundings, thereby falling more easily. In any case, these terms need to be precisely tested in the computations of stellar models.

Let us now discuss the effect of an external turbulence on the thermohaline convection. As already pointed out by Denissenkov & Pinsonneault (2008), the presence of horizontal turbulence decreases the efficiency of the thermohaline mixing and even suppresses it completely as soon as the horizontal turbulent diffusivity is larger than the thermal diffusivity, in other words as soon as the horizontal Péclet number is larger than one. In the small frame of Figure 1, the curves represent the thermohaline coefficients computed by including horizontal turbulent mixing coefficients of 3000 cm² s⁻¹ (dotted curve) and 15,000 cm² s⁻¹ (dot-dashed curve). These turbulent diffusion coefficients are small, but their consequences for thermohaline convection are important. The basic effect is that the thermohaline coefficient drops for values of R₀ decreasing with increasing turbulence. This means that thermohaline convection needs a larger μ-gradient to take place. If the external turbulence is due to rotational mixing as modeled by Zahn (1992) and co-workers, the importance of the horizontal mixing, much larger than the vertical one, cannot be neglected. As discussed by Zahn (1992), whereas the Péclet number is always smaller than one in the vertical direction, it may become high and even larger than one in the horizontal direction. In this case, simply adding the turbulent diffusion coefficients describing the thermohaline convection and the rotation-induced mixing strongly overestimates the global mixing efficiency.

We have derived a new expression for the thermohaline diffusion coefficient in stars, which occurs in the case of mean molecular weight inversion, taking into account the influence of radiative accelerations and external turbulence. We find that the corrections induced by radiative accelerations are small but should not be neglected in the computations of stellar models. When relatively abundant elements like iron are strongly pushed up by radiation, the momentum transfer to the surroundings leads to a negligible effect on the stellar gas as a whole. The only important influence of radiative accelerations on thermohaline convection is due to the difference in the radiatively induced momentum inside the fingers (falling blobs) and outside (surrounding medium). This difference is small compared to the other effects, but in any case the influence of radiative accelerations on the thermohaline efficiency should be systematically checked in all computations of thermohaline convection induced by the radiative levitation of specific elements.

These results are important as they confirm the idea that the accumulation of radiatively levitated elements in stars leads to mixing in a straightforward way. Up to now, atomic diffusion including selective radiative accelerations was known to qualitatively account for a large number of stellar observations (e.g., chemically peculiar stars). However, atomic diffusion alone leads to abundance variations much too important for a quantitative fit to the abundance determinations. Macroscopic motions such as rotation-induced mixing, mass loss, etc., had to be invoked to try to explain the observations. Here we show that a fundamental process was forgotten in all previous computations: thermohaline instabilities directly induced by the atomic diffusion process itself. All the computations of radiative levitation of chemical elements will have to be revised to take this effect into account. Note that in the case where heavy elements are levitated, helium always sinks due to gravitational settling, as the radiative flux is never able to sustain it in its original abundance. This creates a stabilizing μ-gradient which opposes the destabilizing one due to the heavy elements. All these effects have to be treated together, and the expected result is that heavy elements may still accumulate in connection with the helium depletion, but not as much as obtained with pure diffusion (cf. Théado et al. 2009).

Furthermore, these processes have to compete with external turbulence and particularly rotation-induced mixing. This will be somewhat complex, as it depends on the ratio of the horizontal to the vertical mixing coefficients. We have shown that a large horizontal turbulence decreases and may even suppress the efficiency of the thermohaline mixing. On the other hand, vertical turbulence increases this efficiency. A combination of all these effects, including their close interaction, will have to be studied in real stellar cases.

S.V. acknowledges fruitful discussions with Dr. Haili Hu at the Kavli Institute for Theoretical Physics, University of California, Santa Barbara, during the program "Asteroseismology in the Space Age," 2011 September to December.