Chemical element transport in stellar evolution models

Stellar evolution computations provide the foundation of several methods applied to study the evolutionary properties of stars and stellar populations, both Galactic and extragalactic. The accuracy of the results obtained with these techniques is linked to the accuracy of the stellar models, and in this context the correct treatment of the transport of chemical elements is crucial. Unfortunately, in many respects calculations of the evolution of the chemical abundance profiles in stars are still affected by sometimes sizable uncertainties. Here, we review the various mechanisms of element transport included in the current generation of stellar evolution calculations, how they are implemented, the free parameters and uncertainties involved, the impact on the models and the observational constraints.


Introduction
Almost a century ago Eddington wrote [1]: 'It is reasonable to hope that in a not too distant future we shall be competent to understand so simple a thing as a star'.
During this time, the theory of stellar evolution has been developed and established, and its main predictions confirmed by a large number of empirical tests that have involved photometric, spectroscopic and asteroseismic observations (and solar neutrino flux measurements, after the discovery of neutrino oscillations). Results obtained from stellar model computations are nowadays widely used to develop a vast array of techniques to estimate distances, ages, star formation histories and chemical evolution of star clusters and galaxies, both resolved and unresolved [2].
Obtaining this kind of information from techniques rooted in stellar evolution calculations is a crucial step to address problems like understanding the mechanisms that drive the formation and evolution of galaxies. The accuracy of results gathered from stellar population analyses is tied to the accuracy of the current generation of stellar models; in this respect, one particularly  thorny issue is how to treat the transport of chemical elements in stellar evolution calculations. The problem is that mixing and element transport do not arise from the solution of the equations of stellar structure and evolution, instead they have to be 'added' following recipes that often-as we will seeinvolve a number of free parameters and/or are subject to sizable uncertainties.
On the other hand, the temporal evolution of the chemical abundance profiles within stellar models is a main evolutionary driver, and can be in principle tested through spectroscopic observations of photospheric abundances of key elements, asteroseismic observations (study of non-radial pulsations of stars, which can test thermodynamical properties of stellar interiors that are also affected by the chemical composition), and more indirectly through the effect on star counts (sensitive to evolutionary timescales) and evolutionary paths in colour-magnitude-diagrams (CMDs) or Hertzsprung-Russell diagrams (HRDs), which are all affected by the internal chemical profiles.
The aim of this review is to discuss the various mechanisms of chemical element transport included in the current generation of stellar models, their effect on the evolutionary properties of the models, and the various prescriptions found in the literature, that often produce very different results. This will allow the reader to appreciate the main uncertainties involved, and what properties of stellar models are affected the most.
We start in §2 with a brief summary of the equations of stellar structure and evolution, and an overview of the main evolutionary properties of models of different initial masses, to set the stage for the discussion that follows. The next sections discuss the implementation in stellar modelling, the associated uncertanties and the observational constraints of processes like convection ( §3), semiconvection ( §4), thermohaline mixing ( §5), atomic diffusion ( §6), phase separation upon crystallization ( §7) and rotationally induced element transport ( §8). Section 9 discusses briefly an example of how a combination of several of the mechanisms discussed in the previous sections can explain the puzzling trend of photospheric Li abundances with effective temperature in open clusters, and conclusions follow in §10.

Stellar model computation
With stellar modelling we mean here the calculation of the run of physical (i.e. luminosity, temperature, density and specific heats) and chemical quantities, from the centre to the photosphere of a star of a given initial mass and chemical composition, and their evolution with time. Despite tremendous advances in computing power and computational techniques during the last decades, a full detailed modelling of a star by solving the equations of radiation hydrodynamics in three dimensions (3D) is still unfeasible, and this will be the case for the foreseeable future. This inability to model stars with multidimensional radiation hydrodynamics is a consequence of the extreme range of spatial and temporal scales 1 that need to be resolved when calculating full evolutionary models covering all stages from the pre-main sequence (pre-MS) to the last white dwarf (WD) or pre-supernova phases. Current hydrodynamics computations are, however, starting to be able to provide some guidelines about mixing processes [3], as highlighted in the following sections.
For these reasons, complete stellar evolution computations still have to rely basically on the following 'classical' set of one-dimensional (1D) equations for spherical, non-rotating and non-magnetic stars (equation of continuity of mass, hydrostatic equilibrium, 2 energy generation and energy transport, respectively) and Raphson-Newton solution methods [4]: where the independent variable m is the mass enclosed within radius r, and T, L, P and ρ are, respectively, the temperature, luminosity, pressure and density at the layer specified by the value of m. The coefficient ν denotes the energy per unit time and unit mass carried away by neutrinos (that do not interact with the stellar gas), n is the energy per unit time and unit mass produced by nuclear reactions. For a generic nuclear reaction n A A + n b b →products involving elements A and b with mass fractions X A and X b and atomic weights A A and A b where Q Ab is the amount of energy released by a single reaction and R Ab is the number of reactions per unit mass and unit time, given by (2.5) where σ v Ab is the reaction cross section. 3 The coefficient g represents the so-called gravitational energy produced per unit time and unit mass, and is given by where U is the internal energy per unit mass, v = 1/ρ is the specific volume and μ is the mean molecular weight of the stellar matter. The term (∂U/∂μ) T,v (∂μ/∂t) arises from the variation of U at constant temperature and volume due to the change of chemical abundances. Its contribution to the energy budget is negligible when nuclear reactions are efficient, but is very important in the case of WDs, where nuclear burnings are inefficient. When integrated over the whole stellar structure, g is equal to the time variation of the internal energy plus the gravitational potential energy of the star [2,5].
For the case of radiative plus electron conduction energy transport, the gradient ∇ ≡ (d ln(T)/d ln(P)) is set to ∇ rad where a is the radiation density constant, c is the speed of light, G is the gravitational constant and κ is the Rosseland opacity, including also the contribution of electron conduction when appropriate. In the case of convective energy transport, a theory of convection is needed to calculate the appropriate temperature gradient ∇ conv (see §3). These equations are complemented by a set of I equations (s = 1, . . . , I) for the change of the mass fraction of the I chemical elements considered at the layer specified by m. Consider first the changes due just to nuclear reactions; an element s is produced by w reactions of the type n h h + n k k → n p s and destroyed by l reactions n d s + n j j → n z z. (2.8) The complete system of equations (equations (2.1)-(2.4) and I-times equation (2.8)) is solved at a given time t considering the 'Lagrangian' independent variable m, with r, L, P, T as unknowns, once the stellar mass and initial chemical composition are specified, and prescriptions for the equation of state of the stellar gas, Rosseland mean opacities, nuclear reaction cross section and energy generation rates, and neutrino production rates are provided. The chemical composition is usually denoted by X, Y, Z, that correspond to the mass fractions of H, He and all other elements collectively called 'metals', respectively. A relative distribution of the metal abundances also needs to be specified. The exact values of the mass ranges depend on the initial chemical composition and the details of the adopted element transport mechanisms. We show also some final product of the evolution of interacting binaries (Novae, Type Ia supernovae). He-core WDs from single stars can form only on timescales much longer than a Hubble time, but are produced nowadays by interacting binary systems. There are observational indications that massive stars with initial mass above approximately 20M might not explode as supernovae, rather collapse directly to a black hole [6].
Notice that the chemical abundance profile enters explicitly equation (2.8) and n , and affects the coefficients ν , g , the opacities κ and the equation of state, which all depend on the chemical composition of the stellar matter. Figure 1 shows an overview of the evolutionary paths of single stars with different initial masses (plus the main by-products of interacting binary evolution), as derived from complete stellar evolution models. This general evolutionary framework is solid and does not depend on the details of element transport modelling although the precise values of the various mass ranges do (and they depend also on the initial chemical composition). It constitutes a reference guideline for the discussions presented in the following sections.
It is clear that chemical element transport does not arise naturally from the equations of stellar structure and evolution. This makes it necessary to first identify all possible mechanisms-in addition to the nuclear reactions-able to change the chemical abundances at a given mass layer, and then to develop formalisms for their implementation in the 1D equations (a major difficulty). The element transport mechanisms described below will add extra terms to the right-hand side of equation (2.8), in addition to the terms describing the effect of nuclear reactions.

Convection
Besides radiation and electron conduction, convection is the third fundamental mechanism for energy transport in stars. It involves organized large-scale motions of matter that in addition to carrying energy are also a very efficient source of mixing. It can be envisaged as a flux of matter from deeper-hence hotter-stellar layers moving vertically outwards into cooler layers, and material from cooler outer layers flowing down to hotter inner layers. The implementation of this element (and energy) transport in stellar models requires first a criterion for the onset of the convective instability, and then a mathematical treatment to predict the main physical properties of convective regions.

Instabilities in non-rotating stars
Matter inside the stars is never at rest, but usually the gas is locally subject to small random perturbations around equilibrium positions. Under certain conditions, these small random perturbations can trigger large-scale motions that involve sizable fractions of the total stellar mass. These large-scale motions are called convection, and are the equivalent to the motion of water elements in a kettle heated from below.
The treatment of convection in stellar interiors is extremely complicated and requires the introduction of various approximations. This stems from the fact that the flow of gas in a stellar convective region is highly turbulent, forcing us to adopt simplified models that provide only mean approximate values for the properties of the flow of gas. The so-called mixing length theory (MLT) is the local, time-independent convection model almost universally used in stellar evolution calculations, which we will discuss below. Firstly, we show how a simplified linear analysis is sufficient to determine the main criteria for the onset of mixing (not only convective) in stellar interiors [7].
We consider a gas element at rest at a distance r from the star centre. This gas bubble will have a pressure P 0 , temperature T 0 , density ρ 0 and mean molecular weight μ 0 (the molecular weight is the mean mass of the gas particles in atomic mass units) equal to those of the environment, supposed to be in radiative equilibrium (in this section 'radiative' actually means 'radiative plus conductive') as shown in figure 2. If random motions displace the bubble by a small amount r away from the equilibrium position, the equation of motion for an element of unit volume can be written as (assuming the viscosity is negligible) where ρ is the difference (ρ bubble − ρ surr ) between the bubble (supposed to have constant density) and the surroundings, and g is the local acceleration of gravity. One reasonable assumption is that the motion of the bubble is fast enough that all time derivatives of the mean stellar properties are equal to zero. As a consequence of the displacement from its equilibrium position, two distinct physical situations can occur, as shown in figure 2. If the region is convectively stable, the displaced gas parcel experiences a restoring force that moves it back towards the original position, as illustrated in figure 2b. Indeed, the blob is subject to a stable oscillation around an equilibrium position with a frequency named the Brunt-Väisälä frequency (see below). If the bubble of material at position r + r is less dense than the surrounding material, it will continue to be pushed upwards by buoyancy forces and the region is then said to be convectively unstable, as shown in figure 2c. Similarly, if the bubble displaced at a radial distance r − r is denser than the surrounding material, its motion will continue because it is heavier than the local environment. To determine which solution applies to a given layer within a star, we assume that along the displacement r the bubble is always in pressure equilibrium with the surroundings, that is, P = (P bubble − P surr ) = 0, and that the molecular weight of the bubble μ bubble is always equal to its initial value μ 0 (there is no matter exchange with the surroundings, hence the gas parcel retains its identity).
The assumption of pressure equilibrium means that the motion of the bubble has to happen with a speed lower than the local sound speed. When there is a molecular weight gradient dμ/dρ throughout the region, the difference μ = (μ bubble − μ surr ) will be equal to μ = μ 0 − [μ 0 + (dμ/dr) r], hence μ = −(dμ/dr) r.
Using the relationships d ln(P) = (1/P) dP and d ln(μ) = (1/μ) dμ, one gets Differentiating with respect to time, one obtains The temperature difference T = (T bubble − T surr ) depends on the difference between the temperature gradients in the bubble and in the surroundings and the rate of temperature change due to energy losses from the bubble (due, for example, to thermal diffusion) whose efficiency will depend on a parameter we denote as ζ , hence By introducing the notation ∇ ≡ (d ln(T)/d ln(P)) and differentiating with respect to time, one obtains If P = 0, and in the assumption that the differences T, ρ and μ are small, we obtain from the equation of state We have derived in this way the following set of four homogeneous equations for the four unknowns T, ρ, μ and r: and One can search for solutions of the form x = Ae nt . By inserting into the respective equations this functional dependence for T, ρ, μ and r, a non-trivial solution is found when the determinant derived from the coefficients of A T , A ρ , A μ and A r is equal to zero, giving The condition that at least one of the n is real and positive or imaginary with a positive real part (unstable solutions) is given by the Hurwitz criterion, resulting in at least one of the following conditions  to be satisfied (we recall that the pressure P always increases towards the star centre, χ μ is negative, χ T and χ ρ positive): and ∇ rad > ∇ ad . (3.8) If μ < 0 (μ increasing towards the surface), the medium is always unstable. When the gas parcel is displaced upwards (downwards) by a small distance r, its density will be lower (higher) than the environment, and will continue to be pushed upwards (downwards) by buoyancy. The temperature difference between the displaced mass element and its surroundings suppresses or favours this displacement, depending upon the difference between ∇ rad and ∇ ad (more below). For the more common case of μ ≥ 0 (heavier elements are usually synthesized by nuclear reactions in the central regions), we should consider equation (3.7)-the so-called 'Ledoux criterion'-and equation (3.8)-the so-called 'Schwarzschild criterion'. The second term on the right-hand side of equation (3.7) is positive for a positive ∇ μ (the composition gradient has a stabilizing effect) so that if the gradients satisfy the Ledoux criterion, the Schwarzschild criterion is automatically satisfied, hence the Schwarzschild criterion determines the presence of an unstable medium.
These instability criteria are 'local', in the sense that they can be applied layer-by-layer without accounting for non-local effects that can be, however, relevant when dealing with convective mixing. Figure 3 displays a qualitative sketch in the ∇ μ − (∇ rad − ∇ ad ) diagram of the region where instabilities occur, divided into four quadrants. The only stable region in this diagram is the bottom left quadrant, where ∇ μ ≥ 0 and ∇ rad < ∇ ad . All other quadrants are unstable regions, although the type of mixing involved depends on the exact values of the gradients. If ∇ μ ≥ 0 and ∇ ad < ∇ rad < ∇ L (lower right quadrant), the instability is called 'semiconvection'.
The linear analysis in the semiconvective regime shows that if the gas bubble is displaced upwards and ζ =0, its internal temperature will be slightly larger compared to the surroundings, hence its lower density will favour a continuation of the displacement. On the other hand, μ of the raising bubble will be larger than the environment and in the semiconvective regime of the gradients this effect will prevail, causing the return of the mass element to its starting position. When ζ > 0, the temperature of the element returning from above is smaller than its initial value at r, therefore the restoring force is larger upon returning than when the bubble was leaving its starting position. As a consequence, the bubble returns back to a radial location r with a larger velocity than when it left for the upwards motion. On the downward excursion, the same effect (with an opposite sign) will operate and the amplitude of these oscillations will progressively increase due to the increasing radiative losses during the oscillatory motion.
The growth of these oscillations (overstability) will lead to chemical mixing and thus decrease or destroy the stabilizing gradient ∇ μ . Results from numerical simulations and laboratory experiments show a more complex picture [8,9], with in some cases, the formation of well-mixed fully convective layers, separated by stratified interfaces, but the linear analysis suffices to highlight the peculiarity of the mixing associated with semiconvection.
The efficiency (timescale) of semiconvective mixing is difficult to estimate and depends on the efficiency of the bubble energy losses; an analysis of the growth rate in the framework of the linear analysis shows that as a result of semiconvective mixing ∇ rad ∼ ∇ ad [10], and this is another approach to treat semiconvective mixing in stellar evolution calculations, i.e. to impose a mixing efficiency such that ∇ rad ∼ ∇ ad (actually this was the original way in which semiconvection was implemented [11]) in a semiconvective region (see §4).
If ∇ μ < 0 and ∇ rad < ∇ L (top left quadrant), the instability is usually called a 'thermohaline' instability. If ζ = 0, a gas bubble displaced downwards will have a larger μ than the environment, but also a larger temperature than the radiative environment. The combined effect is that the bubble density is lower than that of the environment, and buoyancy will push back the displaced gas element. However, if ζ > 0, the energy loss will eventually decrease the bubble temperature enough to induce a further displacement downwards due to the effect of the larger μ.
This instability is controlled by the heat leakage of the displaced element, and it is observed, for example, when a layer of hot salt water lies over a layer of cold fresh water. Upon displacing downwards a blob of hot salt water, due to the fact that heat transfer by molecular collisions is typically faster than the motion of chlorine and sodium ions that cause the composition to equilibrate, the sinking blob will be able to come into thermal equilibrium with the surrounding medium faster than achieving composition equilibrium. Given that salt water is heavier than fresh water at the same temperature, the blob will continue to sink in the surrounding fresh water. As this motion continues, the medium develops 'fingers' of salt water reaching down into the fresh water. Inside a star, the role of salt can be played by a heavier element like helium, in a hydrogen-rich medium. This is an example of so-called 'doubly diffusive instabilities', because it involves the diffusion of two different components (particles and heat).
In general, a convectively unstable region mixes matter on very short timescales compared with evolutionary timescales (either nuclear or Kelvin-Helmholtz timescales), and the chemical profile in convective layers can always be assumed uniform to a very good approximation (instantaneous mixing approximation), with abundances of individual elements equal to values averaged over the whole convective region. If the convective region extends from mass layer m 1 (inner boundary) to mass layer m 2 (outer boundary) within the star, inside this region the abundance X s of a generic element s is constant. At the boundaries (one or both of them), one may have a discontinuity between the homogeneous convective chemical profile and the profile in the radiative regions, for example, due to nuclear reactions (previous and/or present). Owing to these effects and just to give an example, the time evolution of X s (in the approximation of instantaneous mixing) within an expanding convective shell is to a first order given by where m = m 2 − m 1 , and X s1 , X s2 are the abundances on the radiative side of the discontinuities at, respectively, the inner and outer boundary of the convective region. The first term in the integral describes the variation due to the nuclear burnings (if efficient), whereas the other two terms describe the change in composition when the boundaries of the convective zone move into surrounding regions of-in principle-inhomogeneous composition.
Exceptions to the validity of the instantaneous mixing in convective regions are the advanced evolutionary phases of massive stars about to explode as Type II supernovae (figure 1), the production of Li due to the Cameron-Fowler mechanism in the envelopes of red giant branch (RGB) stars [12] and proton ingestion episodes into the intershell convection zone of low-metallicity asymptotic giant branch (AGB) stars [13,14]. In this case, the nuclear-burning timescales in the convective regions are comparable to convective mixing timescales. This case is usually treated with a time-dependent convective mixing discussed in the next subsection.
Before closing this section, we just introduce a quantity that will appear often in the rest of the paper. In the case of ζ = 0, gas elements in a stable region will oscillate around their equilibrium position with . Schematic illustration of the MLT approximation to convective motion. The mixing length Λ corresponds to the characteristic radial distance scale over which rising and falling convective elements move before merging with the surrounding medium.
a frequency called the Brunt-Väisälä frequency, usually denoted with N. This can be derived easily from equation (3.5), and is equal to Another expression for this frequency (that we will use in the rest of the paper) involves the derivatives δ = −(∂ ln ρ/∂ ln T) P,μ and φ = (∂ ln ρ/∂ ln μ) P,T , and the pressure scale height H P ≡ −dr/d ln(P) = P/(ρg) (where g is the gravitational acceleration), and is Using δ and φ, the gradient ∇ L (Ledoux gradient) can be rewritten as

The mixing length theory of convection
From the point of view of the evolution of chemical abundances, the instantaneous mixing approximation does not require a model for stellar convection. However, this is necessary for calculating convective velocities and fluxes, and when a time-dependent description of convective mixing is required. The formalism almost universally used in stellar evolution calculations is the MLT [15], a simple, local, time-independent model, firstly applied to stellar modelling by Biermann [16]. The formulation by Böhm-Vitense [17] is usually employed in modern stellar evolution calculations. The basic idea of the MLT is to assume that the stellar fluid is composed of identifiable convective elements that move vertically in the gravitational field between regions of higher and lower temperature. Indeed, there is no net mass flow, but the effect is an outward transport of energy. The MLT assumes a characteristic distance over which bubbles rise before dissipating, the so-called mixing length Λ (figure 4), and describes the motion of these bubbles over the characteristic scale Λ under some general assumptions: -all bubbles have the same characteristic size that is of the same order as Λ; -Λ is much smaller than any other length scale of physical significance in the star; -the physical properties, i.e. temperature, density, pressure and chemical composition, of the bubbles differ only slightly from the surrounding medium; -pressure equilibrium with the environment is maintained. This means that the velocities of the convective elements are small compared with the local sound speed in the local environment. Λ is assumed to be equal to a multiple of the local pressure scale height H p , i.e. Λ = α MLT H P , with α MLT -the mixing length parameter-being a constant to be empirically calibrated (usually by reproducing the solar effective temperature at the solar age with a solar model).
The need for a mean free path in the simple framework of the MLT can be easily explained as follows. Let us consider the cross sections of the rising and falling gas columns; if originally in a given convective layer the cross sections were the same, the rising gas (always in pressure equilibrium with the surroundings) will expand by a factor e after a distance equal to H P . This means that at this point within the star there is much less space available for the falling gas. On the other hand, the amount of falling material must be the same as the rising one, otherwise the star would either dissolve or concentrate all its mass in the interior, thus violating the hydrostatic equilibrium condition. The only solution is that after a distance Λ of the order of H P , part of the material stops and inverts its motion.
By relying on the previous assumptions, the MLT provides the following equations for the velocity of the convective elements v c , the convective flux F c and the convective efficiency Γ (defined as the ratio between the excess heat content of a raising convective bubble just before its dissolution, and the energy radiated during its lifetime) [18]: and where ∇ is the temperature gradient of a rising (or falling) element of matter within the convective region, and ∇ c is the average temperature gradient of all the matter at a given level within the convective zone (the quantity needed to solve the stellar structure equations). There are three additional free parameters besides α MLT , i.e. a, b and c, that are usually fixed a priori and define the MLT 'flavour' [19] 5 There are four unknowns, namely v c , F c , ∇ c , ∇ and the three previous equations plus an additional equation arising from the fact that the total flux to be transported by radiation plus convection is known from the solution of the stellar evolution equation If convection is efficient in the deep stellar interiors, the MLT provides ∇ → ∇ ad and velocities of the order of 1-100 m s −1 , many orders of magnitude smaller than the local sound speed. On the contrary, in convective layers close to the surface the gradient is strongly superadiabatic and velocities are much larger, of the order of 1-10 km s −1 , close to the local sound speed.
External or inner regions of a stellar model are convective in the following cases: -Large values of the opacity κ. Given that ∇ rad ∝ κF, where F is the energy flux, ∇ rad tends to increase above ∇ ad . The radiative opacity generally increases with decreasing temperature-for a given chemical composition-hence this situation occurs most commonly in the cooler outer layers of stars. -If the energy generation rate in the star is very sensitive to the temperature, then the energy flux F rises rapidly as r approaches zero in the stellar centre. This large heat flow can eventually cause ∇ rad to increase above ∇ ad . This situation occurs only in stellar cores, when the nuclear energy generation rate is very sensitive to temperature, as is the case of the H-burning CNO-cycle, Heburning and more advanced nuclear-burning stages. -In ionization zones the adiabatic gradient can decrease below the typical value ∇ ad ∼ 0.4. Also, in these regions the radiative opacity usually increases. Owing to these effects, one can expect the ionization zones located in the outer layers to be convective.
When the convective mixing timescales are comparable to nuclear-burning timescales, a diffusive approach is usually employed to follow the chemical evolution [20]. This means that an extra term  is added to the right-hand side of equation (2.8) for a given element s, i.e. the right-hand side of the following diffusion equation: where the diffusion coefficient associated with the convective transport is taken to be D c = 1 3 α MLT v c H P , using the value of v c derived from the MLT (see [21][22][23], for examples).
The MLT is very appealing for its mathematical simplicity (despite the free parameters involved) and the fact that it relies just on local quantities, hence it is easy to include in stellar evolution codes and does not affect the stability of the numerical solution of the stellar evolution equations. Once α MLT is calibrated on the Sun, the resulting HRDs (and CMDs) of the models generally reproduce the observations satisfactorily within the current errors. 6 Comparisons of the effective temperature evolution of models calculated with a solar α MLT and with a calibration of α MLT obtained from 3D radiation hydrodynamics simulations of stellar envelopes [25] display an agreement within 30-50 K [26]. 7 On the other hand, helioseismic data (the study of non-radial oscillations of the Sun [28]) clearly point to shortcomings of the MLT description of the physical structure of convective regions [29].
An alternative local description of convection included in stellar evolution calculations with the ATON code [30] is presented in [31]. The main difference from the MLT is that the convective elements have a spectrum of sizes, and the scale length of the convective motions is set to be equal to the distance to the closest convective boundary. Stellar evolution tracks calculated with this convection model show a different evolution of the effective temperature T eff compared to calculations with the MLT, especially along the RGB phase. Differences increase with increasing initial metallicity. This convection model also suffers from shortcomings when compared with helioseismic results [29].
There are very sophisticated non-local Reynolds stress models that describe not only convection, but also (see below) semiconvection and overshooting [32][33][34][35][36][37]. They introduce a large number of equations to be coupled to the stellar structure equations and several free parameters to be calibrated observationally, and are generally not included in current stellar evolution modelling.
Convection plays a major role in determining the evolutionary properties of stellar models. Figure 5 shows the HRD of models with (the 5M model) and without (the 1M model) convective cores during   In the inner radiative layers of the 1M model, the H-abundance profile changes smoothly during the MS evolution, and at any given time the H abundance increases gradually from the core towards the more external layers where the burning becomes progressively less efficient. In the convective inner layers of the 5M model-where the burning takes place-the H-abundance profile is uniform, and the progressive retreat with time of the outer border of the convective core produces at a given time t a H-abundance profile flat in the innermost regions, then increasing outwards. Surface convection (present everywhere along the HRD evolution, apart from the MS phase of the 5M model) is very important when comparing the surface abundances measured from spectroscopy with the models, because of the dredge-up phenomenon, whereby the fully mixed convective envelopes reach layers processed by nuclear burnings, hence altering their chemical composition during the RGB and AGB phases [2].

Convective overshooting
As already mentioned, the criteria for the onset of convection are local. The boundaries of the convective regions are fixed at the layer where the random motions of the gas do not get amplified. At this convective border, the acceleration of the gas elements is zero, but not the velocity. One expects, therefore, the chemically mixed region to extend into the formally stable layers in the case of both core and envelope convection. In real stars mixing beyond the formal boundary most probably results from the interplay of several physical processes [38][39][40], grouped in stellar evolution modelling under the 'umbrella' term overshooting (or sometimes convective boundary mixing). These overshooting processes are modelled very crudely by introducing free parameters to be calibrated on observations of eclipsing binaries [41,42] (typically comparisons with masses and radii of both components under the assumptions that they are coeval) and open clusters [43] (the shape of the MS turn-off (TO) region).
The standard approach is to assume that overshooting into the stable layers does not affect their thermal structure, hence the temperature gradient stays radiative. The composition is then mixed instantaneously between the formal convective border and layers at a distance λH P from this border, where λ (λ < 1) is a free parameter to be observationally calibrated, and H P is the pressure scale height at the convective border. This is the approach taken, for example, in the BaSTI [44], DSEP [45], Yale-Yonsei [46] codes.
A similar approach of instantaneous mixing is taken in [47]. The following integral criterion is employed to fix the extent of the convective core overshooting, inspired by a constraint on the maximum possible extension of the overshooting region developed by Roxburgh [48] where F over is a free parameter (between 0 and 1) to be calibrated against observations. The two luminosities L rad = −((16π acT 3 r 2 )/(3κρ)) dT/dr and L N represent, respectively, the local radiative luminosity and the total luminosity produced by nuclear reactions within radius r, while r cc is the radius of the convective core boundary and r ov is the radius of the outer boundary of the overshooting region, to be determined through this equation. The various releases of Padua stellar evolution models [49] and the PARSEC code [50] consider instead an adiabatic gradient for the overshooting region (this is the case of penetrative convection according to [51], because the overshooting material is assumed to be able to change the entropy stratification), and the spatial size of the overshooting region is determined as follows [52]. Starting from each radial distance r i from the centre inside the formally convective region, the following equation (based on the MLT) is integrated outwards, up to r = r i + l, where l = λH p , H p being the pressure scale height at the convective boundary, and λ (λ < 1) a free parameter: is the total energy flux and F rad = −((4acT 3 )/(3κρ)) dT/dr. In the overshooting region, where the actual gradient is set to adiabatic, F rad > F hence F c < 0, mimicking the fact that the convective elements penetrating into the stable layers are formally cooler than the surrounding medium. At any given r, the convective elements originated in the range between r − l and r will display a range of velocities, and the maximum value v m attained at each r is taken to derive the run of v m (r) as a function of r. The border of the instantaneously mixed region is then taken at the radial distance for which v m = 0.
In all these approaches, the free parameters that determine the spatial extent of the convective core overshooting region needs to decrease in the regime of small convective cores. The reason is that H P → ∞ as r → 0, hence the smaller the size of the formally convective core, the larger the size of the actual mixed core (including overshooting region), with no smooth transition of the shape of MS evolutionary tracks from the convective core to the radiative core case. Observations of the MS TO of open clusters of various ages confirms the need to decrease to zero the extent of the overshooting regions when the stellar mass decreases, in the regime of small convective cores (masses below approx. 1.5M -see, e.g. [47]). This means that the assumed trend of λ (or F over ) with mass, for masses with small convective cores, introduces an additional degree of freedom, albeit affecting only a restricted range of stellar masses. In this context, following an integral constraint on the maximum possible extent of overshooting regions [48,53], it has been proposed to limit the extent of the overshooting region to 15% of the radius of the formally convective core to generate a smooth transition from models with convective cores to models with radiative cores on the MS [54].
Additionally, one can find in the literature two very different approaches to the mixing beyond the formal convective boundary.
(i) A diffusive approach [55] used, for example, in the STARS [56], MESA [57,58] and GARSTEC [59] codes. The element transport beyond the formal convective boundary is described as a diffusive process (avoiding the instantaneous mixing approximation) based on results of 2D radiation hydrodynamics simulations of shallow stellar surface convection zones [60] The diffusion coefficient D ov is given by where D c is the diffusion coefficient inside the convective region (D c = ( 1 3 )α MLT v c H P ), z is the distance from the convective boundary, H P is the pressure scale height at the convective boundary and f is a dimensionless free parameter. Typical calibrated values of f are approximately 0.01-0.02. (ii) Modelling of the mixing as 'turbulent entrainment', following the simulations by Meakin & Arnett [61], as employed in the stellar evolution calculations by Staritsin [62]. The motion of matter in the zone of convective instability is turbulent, and the rising turbulent flow spreads horizontally near the boundary of the turbulent region. The interface between the convective region and the stable layers moves through the stable region due to the continuous involvement of new layers in the turbulent motion. The velocity V e of the penetration of the convective turbulent boundary into the stable layers is determined by the 'turbulent-entrainment' law written as where V t is a typical turbulent velocity at the boundary (that can be taken, for example, from the MLT), A and n are parameters that characterize the entrainment (A ∼ 0.027 and n ∼ 1.05 according to the simulations by Meakin & Arnett [61]), and Ri B is the so-called bulk Richardson number defined as with l denoting the typical size of the eddies doing the entrainment (some fraction of the pressure scale height H p at the boundary of the mixed region) and This integral is performed across a region of thickness h that contains the convective boundary, N 2 being the the Brunt-Väisälä frequency. When V e is determined, the distance d over which this boundary shifts during an evolutionary time step t is given by d = V e t.
Despite the uncertainties involved in its parametrization, the treatment of overshooting is very important because it affects the evolutionary properties of the models. For example core overshooting during the MS produces brighter models, an increased MS lifetime (because more fuel is available), larger He-core masses that induce a brighter and shorter-lived He-burning phase, and also less extended loops in the HRD (figure 5). Overshooting below convective envelopes can alter the surface abundances after the dredge-up episodes [63], and affect the luminosity of the RGB bump (see §5) and also the extension of the loops in the HRD (in the case of the 5M model in figure 5, the loop during the core He-burning phase would become more extended with overshooting also from the convective envelope). Age estimates of young-intermediate age clusters are obviously affected by the amount of overshooting included in the models, i.e. the larger the overshooting region the older is the age estimate of a given cluster. Overshooting also plays an important role during C-burning in super AGB stars, affecting the propagation of the carbon-burning flame [64]. Too large overshooting at the base of the convective Cburning region can prevent the flame from reaching the centre, thus producing a hybrid CONe core, and eventually a CONe WD.
It has also been shown that the diffusive approach to mixing in the overshooting region provides different results in terms of evolutionary times compared to instantaneous mixing, because of a slower addition of extra fuel from the overshooting region when this scheme is implemented [23].

Semiconvection
After convection, semiconvection (called 'double-diffusive convection' in oceanography) is probably the most significant element transport mechanism in non-rotating stellar evolution models. Below are the two major cases where semiconvective transport is efficient, and how this is implemented in stellar models. 8

H-burning phase with convective cores
There is a large body of literature that addresses the issue of semiconvection in massive stars with convective cores during the MS [10,11,[66][67][68]. In a nutshell, layers left behind by shrinking convective cores during the MS will have a hydrogen abundance that increases with increasing radius (the chemical composition of each layer is determined by the composition of the convective core at the moment the layer has detached from the retreating mixed region). They are characterized by ∇ ad < ∇ rad < ∇ L , and a treatment of semiconvective mixing is needed ( figure 7). There is also a narrow range of masses around 1.5M (the values depending on the initial chemical composition) with increasing mass of the convective core during part of the MS, where a narrow semiconvective region forms right above the fully mixed core [69]. The efficiency of element transport in this semiconvective region has important consequences for the morphology of the TO of the tracks-very efficient semiconvection mimics the case of overshooting beyond the formal convective boundary and no semiconvection [69]-and the post-MS evolution of massive stars [70,71] because it changes the H-profile above the H-exhausted region. To provide a guideline, for massive stars inefficient mixing favours core He-ignition on the cool (red) side of the HRD, while efficient mixing favours the ignition on the hot (blue) side of the HRD (and increases MS evolutionary timescales). One could think that comparisons with the observed ratio of blue to red supergiants (B/R ratio) could provide strong constraints on the efficiency of semiconvective mixing in massive stars. However, other factors like the extent of the overshooting region (that acts in the direction to reduce the size of the semiconvective layers), and rotation will affect the HRD location (and eventually loops) during the He-burning phase [72]. At the moment no set of theoretical models seems to be able to match observational constraints on the B/R ratio, that is also strongly affected by mass transfer in binaries, which make up most of the massive star population.
Semiconvective mixing is included in these evolutionary calculations following various recipes. A traditional approach was to iterate the composition in the individual semiconvective layers until ∇ rad = ∇ ad locally [11,70,71]. More modern calculations adopt again a diffusive approach, including a chosen semiconvective diffusion coefficient D SC . The coefficient from [10] is implemented in the codes STERN [73], MESA and GARSTEC. It is derived from a linear local stability analysis [7], assuming the MLT to determine the velocity of the gas elements (using α MLT = 1.5), and is given by is the thermal conductivity, ∇ is the actual temperature gradient and the free parameter α SC determines the mixing timescale (larger α SC correspond to shorter mixing timescales). Values of α SC currently used are of the order of approximately 0.04-0.1, calibrated on empirical constraints [22,72]. The value of ∇ in the semiconvective region is determined from where β is the ratio of the gas pressure to the total pressure (gas plus radiation), L is the total luminosity and L rad is the radiative luminosity, which can be written as Larger values of α SC produce semiconvective temperature gradients ∇ increasingly close to ∇ ad . Another expression for D SC has been derived in the assumption of layering of the semiconvective region, with nearly uniform composition in each layer, separated by thin boundary layers within which the chemical elements are transported by molecular diffusion alone [74]: where K T = K/(ρc P ) is the thermal diffusivity, β is the ratio of the gas pressure to the total pressure (gas plus radiation) and D s is the diffusion coefficient of He due just to the He abundance gradient (see §6). For a chemical composition of essentially two elements (in our case H and He) with atomic numbers Z 1 and Z 2 , masses m 1 and m 2 and number densities n 1 and n 2 , with λ D = (K B T/(4π n e e 2 )) 1/2 (Debye length) and n e the electron density. Typically, D s is very small (smaller by about eight orders of magnitude) compared with K/(c p ρ) (the so-called thermal diffusivity) and the predicted semiconvective transport is very inefficient.
The code KEPLER [75,76] implements a different diffusion coefficient for the semiconvective transport: where D c is the diffusion coefficient the layer would have if the Schwarzschild criterion is used (fully efficient convection, D c = 1 3 α MLT v c H P ) and α SC is a free parameter usually fixed to 0.1 in KEPLER calculations [77].
Finally, results from recent 3D hydrodynamics simulations of layered semiconvective regions [9] have been transposed into a diffusion coefficient D SC implemented in the MESA code. The result is that in stellar conditions the mixing obtained with this coefficient is very fast and is essentially equivalent to calculations performed with instantaneous mixing in the semiconvective region [78].

Core He-burning phase in low-intermediate-mass stars
Another important evolutionary stage where semiconvection plays a major role is the core He-burning phase of low-and intermediate-mass stars, which has been widely discussed in the literature (see [79][80][81][82] and references therein). If we consider a low-mass star (initial mass below approx. 2M ) after the core He-flash, He-burning is efficient in a convective core with central values log(T c ) ∼ 8.07 ± 0.01, log(ρ c ) ∼ 4.25 ± 0.05, and a chemical composition of almost pure He. The opacity is dominated by electron scattering, but an important contribution (about 25% of the total) comes also from free-free absorption. The transformation of He to C due to nuclear burning increases the free-free opacity, hence κ, and within the convective core ∇ rad increases, developing a discontinuity of ∇ rad at the inner convective boundary ( figure 8).
This discontinuity of the radiative gradient is clearly unphysical. In fact, in convective regions far from the surface ∇ c ∼ ∇ ad , and given that in the MLT F c ∼ (∇ rad − ∇ c ), there would be an increasing convective flux developing on the inner side of the convective boundary. According to [83], the presence of this discontinuity is due to an incorrect application of the Schwarzschild criterion when an MLT picture of convection is considered, for a proper implementation requires always ∇ rad = ∇ ad at the convective side of the boundary. Another interpretation is that even a very small amount of overshooting suffices to reach the neutral ∇ rad = ∇ ad condition at the convective boundary, by means of a self-driving mechanism.
Let us consider an overshooting such that just one radiative layer is mixed. This layer will become fully convective on very short timescales, with locally ∇ rad > ∇ ad because of the amount of carbon mixed from the layers below. This means that it is now the formal boundary of a new enlarged convective core. A small overshooting to mix the next radiative layer will have the same effect, and so on until finally ∇ rad = ∇ ad at the convective side of the boundary. Another possibility is atomic diffusion [84] (see §6), driven by the gradient in carbon abundance between the fully mixed convective core and the surrounding radiative layers. As carbon diffuses into the radiative He-rich layers, the local opacity increases, ∇ rad becomes larger than ∇ ad and these layers become convective.
This diffusion is able to 'extend' the convective core on short timescales (short compared to nuclear timescales), so that ∇ rad = ∇ ad is satisfied at the convective boundary [84]. Yet another possibility to extend the formally convective core to attain ∇ rad = ∇ ad at the convective side of the boundary is the shear instability [85] (see also §8). At the formal boundary of the convective region (where the discontinuity of ∇ rad appears), we have two fluids of mean molecular weight μ 1 and μ 2 , respectively, whose surface of separation is perpendicular to the gravity field, and have a relative velocity v tang tangential to the separation surface. In this situation v tang ∼ v c , the velocity of the overturning convective elements when they reach the formal convective boundary, as obtained from the MLT. 9 In this case, a mixed transition region should appear whose width is By denoting with t mix the typical time to mix this region of width z, and assuming that t mix is of the order of the characteristic time of convection at the boundary of the core (as suggested by Castellani 9 Even a difference in rotational velocity between the convective core and the overlying radiative layers can have a similar effect. In this where v tang is now the difference in rotational velocity between the convective core (expected to rotate like a solid body) and the overlying radiative layers. Equation It has been shown in [79] using an MLT approach that an overshooting length z ≈ 5 cm guarantees an advancement of the fully mixed core with speed v p ≈ 10 3 [1 − (∇ ad /∇ rad )] cm yr −1 . This is sufficient to enlarge the mass of the fully mixed core to the point where ∇ ad = ∇ rad ) on timescales shorter than nuclear timescales. By adopting the values of g, v c , λ and μ 1 /μ 2 adopted by Castellani et al. [79], equation (4.5) and (4.6) provide z ∼ 3 cm, and v p ≈ 10 3 cm yr −1 , consistent with results in [79].
In practical terms, at every computational time step, one can place the boundary of the fully mixed convective region at the layer where ∇ rad = ∇ ad (figure 9).
After this early phase, when the extension of the convective core follows the increase of ∇ rad , the radiative gradient profile starts to show a minimum (profile 3 in figure 10), as a consequence of the progressive outward shift of the convective boundary. The presence of this minimum depends on the complex behaviour of the physical quantities involved in the definition of ∇ rad , such as opacity, pressure, temperature and local energy flux. In this situation, outward mixing to eliminate the ∇ rad discontinuity will induce a general decrease of the radiative gradient in the whole convective core (due to an average increase of He and consequent decrease of C-see profile 4 in figure 10). The radiative gradient will eventually decrease and become equal to ∇ ad at the location of the minimum of ∇ ad (profile 5 in figure 10).
The semiconvective region is the region between the 'neutral' ∇ rad = ∇ rad point and the overlying formally convective shell whose upper boundary still displays a discontinuity in ∇ rad . In fact, a full instantaneous mixing between the convective core located inside the minimum and the external convective shell would have the consequence of decreasing the radiative gradient in the whole mixed region. However, due to the presence of this minimum, ∇ rad in a portion of the mixed core would become lower than the ∇ ad , i.e. it would not be convective. A solution for this inconsistency is to impose a partial mixing in the formally convective shell (see, e.g. [86] for an example of implementation), such that the final chemical composition-shown in figure 11-satisfies the condition ∇ rad = ∇ ad (profile 6 in figure 10). The mass location of the minimum of the radiative gradient moves outwards with time,   because of the evolution of the chemical abundances caused by nuclear burning. As a final result, the C-enriched region increases its size outwards.
The effects of semiconvection on the evolution of these models are the following: the evolutionary tracks perform more extended loops in the HRD, and the central He-burning phase lasts longer because of a larger amount of fuel to burn.
When the central He abundance has decreased to about Y ∼ 0.10, α-captures by C nuclei tend to overcome C production by 3α reactions, thus He-burning becomes mainly a 12 C + α production of oxygen, whose opacity is even larger than that of 12 C. This causes an increase in the size of the semiconvective region and, in turn, more fresh helium is transferred into the core, which is now nearly He-depleted. Even a small amount of He added to the mixed core enhances the rate of energy production, thus the luminosity increases, driving an increase in the radiative gradient. As a consequence, a phase of enlarged mixed zone starts, the so-called breathing pulse. After this breathing pulse, the star readjusts to burn steadily in the core the fresh He driven there by semiconvection. Detailed calculations show that a few breathing pulses are expected before the complete exhaustion of He in the core. The evolutionary effects of the breathing pulses are the following: the models perform a loop in the HRD at each pulse, the He-burning lifetime is slightly increased and the mass of the CO-core at He exhaustion is increased.
Empirical constraints-mainly the number ratio between horizontal branch (HB-core He-burning phase in low-mass stars) and AGB stars in Galactic globular clusters (GGCs)-suggest that the efficiency of the breathing pulses phenomenon is very low, if any [87,88]. Therefore, they are usually inhibited in stellar model computations by using ad hoc numerical assumptions [82,89]. Typically, during the late stages of core He-burning one forces the extension of the mixed region not to lead to an increase in the central He abundance from one model to the next [87]. Another option is to set to zero the gravitational term in the energy generation equation for the inner regions. In this way, the breathing pulses are also effectively inhibited [82].
This brief discussion highlights clearly the difficulty in modelling core mixing in HB stars. Notice that the standard treatment is usually instantaneous mixing, which may not be adequate in a semiconvective regime. None of the diffusive semiconvective formalisms employed for semiconvection related to Hburning is usually applied to this situation, although the diffusive mixing employed in the STARS code still leads to the onset of breathing pulses during the core He-burning phase. It is also worth recalling that with the inclusion of an extended overshooting region (approx. 1 H P ) beyond the layer where the ∇ rad discontinuity develops, the fully mixed core is always so large that the need to include semiconvective layers disappears [90] (breathing pulses still seem to appear also with large overshooting). An important consequence of including semiconvection, or large overshooting, and/or breathing pulses, is that the CO profile in the final CO core changes, with important consequences for the cooling times of the final WD stage. 10 With increasing stellar mass, the weight of the free-free opacity decreases (because of higher core temperatures) and eventually all these problems disappear in the regime of intermediate-mass stars (masses above a few solar masses).

Thermohaline mixing
In recent years, the role played by thermohaline mixing in stellar evolution has been widely explored in connection with low-mass RGB evolution. 11 When the convective envelope deepens after the TO, the surface chemical composition is altered due to the dredge-up of H-burning processed matter-the so-called first dredge-up-that increases the abundance of N and He, and decreases C and the 12 C/ 13 C ratio. The first dredge-up is completed when the convective region reaches its maximum extension, approaching closely (but not reaching) the H-burning shell around the inert (and electron degenerate) He-core. From this moment on the receding convective boundary is not expected to modify further the surface abundances along the RGB. Spectroscopic observations of metal-poor Galactic halo stars provide, however, compelling evidence for an additional mixing process occurring when RGB stars reach the luminosity of the RGB bump (figure 12), which causes a sudden drop of the isotopic ratio 12 C/ 13 C, a decrease in Li and C, and an increase in N. This mixing affects approximately 95% of low-mass stars, At the completion of the dredge-up, the He-core mass is equal to approximately 0.2M (it was approx. 0.14M at the start of the dredge-up) and the H-abundance discontinuity left over by the fully mixed convective envelope at its maximum extension is located 0.25M away from the centre. The varying H abundance in the layers between the He-core and the discontinuity has been produced during core H-burning along the MS. When the H-burning shell reaches the H-abundance discontinuity the RGB model experiences a temporary drop in luminosity (that in a stellar population produces a local increase in star counts, the so-called 'RGB bump' , before evolving again with increasing L when the discontinuity is crossed.
regardless of whether they populate the halo field or clusters [94,95], and thermohaline mixing has been proposed to explain these observations. 12 In low-mass stars, the main H-burning mechanism during the MS is the p-p chain that, due to its weak dependence on temperature, is efficient also in stellar layers far from the star centre. As a consequence, 3 He accumulates in a broad zone outside the main energy production region. During the first dredgeup, this 3 He is mixed within the convective envelope, with the consequence that, during the following RGB evolution, the layers above the H-discontinuity left over by the receding convective envelope at its maximum extension will have a uniform 3 He abundance, larger than the initial one.
When the shell advances towards the surface during the RGB evolution, in the outer wing above the point of maximum burning efficiency there is a narrow region where 3 He is processed through the reaction 3 He( 3 He, 2p) 4 He. In this nuclear reaction two nuclei transform into three and the mean mass per nucleus-the molecular weight μ-decreases. This leads to a small local decrease of μ when moving from the surface towards the centre of the star. As long as the H-burning shell advances through layers below the H-discontinuity-when the star evolves before the RGB bump, see figure 12-this effect is negligible because the shell is moving in a region with a large positive gradient, due to the H profile left over at the end of the MS. However, when the H-burning shell enters the region of uniform H abundance above the discontinuity, the local inversion of the μ profile, of the order of one part in 10 4 , becomes important [97]. This situation corresponds to the conditions for thermohaline mixing (see [98] and references therein).
Thermohaline mixing is usually included in stellar evolution codes as a diffusive process that tends to erase the molecular weight inversion, with a diffusion coefficient derived from a linear analysis [99,100] where C th = ( 8 3 )π 2 α 2 , α being a free parameter related to the aspect ratio (length/width) of the mixing elements. 13 12 About 5% of the post-bump RGB objects investigated do not show these abundance variations. A strong magnetic field can potentially inhibit thermohaline mixing [96]. 13 The calculations by Cantiello & Langer [91] and the MESA calculations by Lattanzio et al. [101] consider the constant C th as C th ≡  As shown in figure 13, evolutionary calculations show that thermohaline mixing extends between the outer wing of the H-burning shell and the inner boundary of the convective envelope, merging with the outer convection in a short time (approx. 30 Myr for a model with mass of the order of 1M ). Therefore, depending on the mixing efficiency (hence the choice of C th ), a significant amount of nuclear-processed matter in the hotter layers of the H-burning shell can be dredged up to the surface during the remaining RGB evolution, helping to explain the spectroscopic data.
Comparisons with observations require C th ∼100−300 from Li observations in GGCs, and C th ∼ 1000 for C in globulars, clearly mutually inconsistent. On the other hand, models in [92,98] are able to match the 12 C/ 13 C isotopic ratio, C and N abundances in halo RGB field stars with C th ∼ 1000; the same choice of C th allows also to match measurements of the carbon isotopic ratio and the [N/C] ratio in open clusters, and Li abundances in field disc stars. Results from 3D hydrodynamics simulations (rescaled to match stellar conditions) add an additional source of uncertainty, for they predict an efficiency equivalent to just C th ≈ 10 [102,103].
A recent detailed numerical analysis [101] has shown that the surface chemical abundances predicted by stellar models accounting for thermohaline mixing depend on numerical assumptions like spatial and time resolution adopted in the computations. As a consequence, the predicted surface chemical abundances should be treated with caution until a firmer assessment on how to treat thermohaline mixing in model computations is achieved.

Atomic diffusion
Microscopic effects related to collisions among the gas particles induce a slow element transport within radiative regions. It is possible to show from first principles that individual ions are forced to move under the influence of pressure as well as temperature gradients, which both tend to displace the heavier elements towards the centre of the star, and of concentration gradients that oppose the above processes. Radiation-which does not have a major effect in the Sun [104]-pushes the ions towards the surface, whenever the radiative acceleration imparted to an individual ion species is larger than the gravitational acceleration. The speed of the diffusive flow depends on the collisions with the surrounding particles, as they share the acquired momentum in a random way. It is the extent of these 'collision' effects that dictates the timescale of element diffusion within the stellar structure, once the physical and chemical profiles are specified.
The most general treatment for the element transport in a multi-component fluid associated with diffusion is provided by the Burgers equations [105]. They are obtained assuming the gas particles have approximate Maxwellian velocity distributions, the temperatures are the same for all particle species, the mean thermal velocities are much larger than the diffusion velocities and magnetic fields are unimportant, and they can be written as dp i dr including the heat flow equations, In addition, there are the constraints of electric current neutrality, In the above 2N + 2 equations p i , ρ i , n i ,Z i and m i denote the partial pressure, mass density, number density, mean charge and mass for species i, respectively. The total number of species (including electrons) is N. The 2N + 2 unknown variables are the N diffusion velocities w i , the N heat fluxes r i , the gravitational acceleration g (the comparison of the derived g with the known value from the integration of the stellar structure equations provides an important check for the consistency of the results) and the electric field E. The coefficients K ij , z ij , z ij and z ij have to be specified, together with the radiative accelerations g rad . Several stellar evolution codes-in one form or another-use the routine by Thoul et al. [106] to solve the Burgers equations and calculate the velocities of the various chemical species. These diffusion velocities can then be inserted as an advection term in the equation for the time evolution of the mass fraction abundance X i where we show on the right-hand side just the contribution of diffusion to the time evolution of the abundance X i of element i. In Burgers' formalism, the effect of collisions between ions is represented by the so-called resistance coefficients, i.e. the matrices K, z, z , z , whose precise evaluation is essential to estimate correctly the diffusion timescales for the various elements. These resistance coefficients can be expressed in terms of the so-called reduced collision integrals Ω (l,s) ij * according to the following relationships: and and 3 (Z i Z j e 2 ) 2 n i n j ln(Λ 2 ij + 1), (6.11) with λ D being the Debye length; Λ ij = 4T * ij the plasma parameter (ln(Λ ij ) is called Coulomb logarithm); μ ij = m i m j /(m i + m j ) the reduced mass; and Z i , m i and n i the charge number, mass and particle number density of species i, respectively.
The collisions between particles of species i, j in the stellar plasma determine the values of the Ω (l,s) ij * integrals; the physics of the collisions is specified by some form of the Coulomb interaction. This, as a first approximation, can be described by a pure Coulomb potential with a long-range cut-off distance, typically equal to λ D . Using this truncated pure Coulomb potential, the resistance coefficients originally computed by Burgers become [105] and where C E is Euler's constant and the ± signs denote repulsive (+) and attractive (−) pairs of particles, respectively. More accurate collision integrals have been obtained by considering a Debye-Hückel type of potential where r is the particle distance. The results are provided either in tabulated form [107] or as fitting formulae [108][109][110][111][112], and some of them are compared in figure 14, as a function of Φ ij = ln(ln(1 + Λ 2 ij )). Given that the calculations shown in the figure (apart from Burgers' results) make all the same physical assumptions, results are similar, but they all differ significantly from Burgers' results when moving towards higher densities (lower values of Φ ij ), where his approximations are no longer adequate.
Some of these authors argued that while λ D is a suitable screening distance at low densities, the Debye sphere loses its significance in denser plasmas, and that a more appropriate screening distance is in this case the mean interionic distance. Hence, they suggested to use the larger value between λ D and the mean interionic distance (in the Sun, the former has always been larger than the latter). Whatever the choice for the actual screening distance, its value has to be employed as λ D in equation (6.10) to compute the appropriate Λ ij for determining the collision integrals. The effect of quantum corrections on the resistance coefficients has been included in [111] while the calculations [112] account for very recent developments in ionic transport properties in strongly coupled plasmas [113].
As for the radiative levitation of element i, the main physical interactions that drive the transfer of momentum from the radiation field to ions are bound-bound and bound-free transitions, while it is practically ineffective for fully ionized elements. The radiative acceleration can be written as [114] g rad,i = μκ μ i c l 4π r 2 γ i (6.17) where μ is the mean atomic weight, μ i is the atomic weight of element i, l is the local luminosity and r is the radius [114]. The dimensionless quantity γ i depends on the monochromatic opacity data  ). Dark black lines denote the case of a repulsive potential, while the brighter ones denote attractive forces (not computed in [108,109]). The dots represent the actual values computed by Muchmore [108], not his fitting formulae.
where u = (hν)/(K B T), σ i is the cross section for absorption or scattering of radiation by element i, a i accounts for the fact that in bound-free transitions only a fraction of the momentum of the ionizing photons is transferred to the ion, the rest being transferred to the electron lost by the ion, and f i is the number fraction of element i. The Opacity Project provides a set of codes called OPSERVER that enables the calculation of g rad,i [115].
Calculations of precise g rad,i values involve carrying out the integration over about 10 4 u values for each atomic species [116]. Given that these calculations have to be repeated at each layer in the stellar model and at each time step during the computation of an evolutionary sequence, this explains why, to date, there are only few extended sets of stellar models that include also the effect of radiative levitation. One has also to note that the Rosseland mean opacity entering the stellar structure equations as well as equation (6.17) has to be continuously recalculated at each mass layer, not just in the nuclear-burning regions-this is in principle true also for calculations including only the other diffusive processes listed above-to be fully consistent with the composition changes.
As an example, figure 15 displays the total radiative acceleration of iron in a stellar envelope. When the radiative acceleration is larger than gravity below a convective envelope-if convection is present, the effect is obviously negligible because of the much faster timescales of convective mixing-the element can diffuse towards the surface. Figure 16 displays the velocity profiles of H (moving upwards), He and CNO (sinking) within a solar model, derived from the solution of the Burgers equations.
The diffusion velocities are always dominated by the effect of pressure gradients. The radiative (upwards) acceleration of C, N and O is only at most about 5% of the local acceleration of gravity just below the convective envelope boundary [104], and decreases fast moving towards the centre. The local sharp increase of the settling velocity of C and the corresponding small decrease for N at r/R ∼ 0.15-0.20 are due to the effect of abundance gradients when these elements attain the equilibrium abundances of the CN cycle (C decreases, while N increases).

The effect of atomic diffusion on stellar models
Atomic diffusion (sometimes denoted as microscopic diffusion, and hereafter simply diffusion) has a major direct effect on the MS evolution, the chemical stratification of the external layers of hot horizontal branch stars 14 and WDs, and the internal chemical stratification of cold WDs. However, some properties of other evolutionary phases are also indirectly affected, as discussed below [84,[117][118][119][120][121].
The impact of diffusion on stellar models will be discussed mainly in the context of low-mass stars, due to their evolutionary timescales comparable to the diffusion timescales during the MS phase. One has also to take into account that massive hot stars, where in principle radiative levitation can be extremely efficient, experience strong mass loss and rotational mixings (see later on in this section and §8) that tend to limit or completely erase the effect of diffusion. If we consider the evolution of a typical low-mass star with a convective envelope on the MS, the surface abundances of metals and He tend to decrease because of diffusion from the bottom boundary of the shrinking convective envelope (that maintains a uniform chemical profile due to the shorter convective timescales compared to diffusion) being replaced by hydrogen. However, if the radiative acceleration on some ion species is larger than the local gravitational acceleration below the convective envelope, these elements are slowly pushed into the convective zone and their surface abundance increases. In general, the variation of the surface abundances during the MS phase is dictated by the interplay between radiative levitation and the sedimentation due mainly to pressure gradients (often denoted as gravitational settling). This variation of the surface abundances reaches a maximum around the TO.
When, after the TO, the convective envelope starts to deepen, it will rehomogenize an increasingly larger fraction of the stellar mass. Once the star reaches the RGB and the first dredge-up is completed, the abundance changes previously developed are almost completely erased, apart for the very inner layers, where metals and He were sinking during the MS. In general, the smaller (in mass) the convective envelope, the larger is the change of surface abundances during the MS, because of the smaller diluting buffer of matter with the initial chemical composition. Figure 17 compares the HRD of low-mass, metal-poor evolutionary tracks, calculated with and without the inclusion of atomic diffusion. The track with diffusion increasingly diverges from the nodiffusion one moving along the MS. Its TO is fainter and cooler, while the two tracks coverge again along the RGB, where the effect of diffusion is practically erased by the fully mixed deep convective zone. The MS evolutionary times are also different, for the track with diffusion has a shorter lifetime because of the slow increase of He at the centre at the expenses of H during the MS phase. 15 These changes of TO luminosity and MS lifetime at a given stellar mass affect also theoretical isochrones, hence age determinations of old stellar populations based on the TO luminosity decrease by about 10%. Figure 18 shows the pattern of variations of TO surface abundances, for a typical stellar model at the TO of a metal-poor globular cluster. These models have thin convective envelopes that maximize the effect of diffusion.
It is straightforward to notice the selective effect of diffusion, due to radiative levitation. Elements like He, Li, C, N and O sink below the convective layers and their surface abundances are severely depleted up to a factor of about ten. Other elements like Mg, Si, S and Fe are pushed into the convective envelope by the radiation pressure and their surface abundances increase during the MS. The size of these abundance variations decrease in models with larger convective envelopes. For example, a 0.8M 15 The surface abundance variations tend to increase the envelope opacities which, in turn, would tend to increase the MS lifetime. But the dominant process is the decrease of available fuel, and the net effect is an approximately 10% decrease of the MS lifetime.  After the first dredge-up is completed, these abundance variations basically disappear. The He abundance post dredge-up is slightly lower in the model with diffusion ( Y ∼0.01), the RGB bump magnitude decreases by approximately 0.07 mag, the He-core mass at the He-flash increases by a few 0.001M . As a consequence, the RGB tip luminosity increases at the level of just approximately 0.01 mag. The abundances of metals (other than CNO) in the He-core are only a few 0.01 dex more abundant than originally, while they (and He) are less abundant than originally by about the same amount around the core.
Moving to the following HB phase, the zero age horizontal branch (ZAHB) luminosity decreases by approximately 0.02 mag compared to the no-diffusion case, driven by the lower amount of He in the envelope, but the major effect is on the surface abundances during the following HB evolution, as shown in figure 19. At T eff above approximately 6000 K, the surface chemical composition of the models starts to be altered by diffusion; the selective effect increases during the HB evolution and is more pronounced for lower HB masses that evolve at higher T eff . The abundance of He tends to decrease, while overall the surface metal content increases. Figure 19 displays, as an example, the predicted behaviour of the surface Fe abundance in metal-poor HB models. One can notice an increase in surface Fe by a factor of approximately 10 within the first 10 Myr, and up to a factor of approximately 100 after 60 Myr.
These surface abundance changes are again erased by the deepening envelope convection for all HB masses that evolve to the following AGB phase. The subsequent hot WD phase sees diffusion altering again in a major way the surface abundances. Owing to the high surface gravity, gravitational settling of metals (down to the edge of the electron degenerate CO core) in the doubly layered H and He (DA type) or He-dominated (non-DA type) envelopes is very fast. Timescales depend on the dominant element, the envelope thickness, T eff and gravity, but they are at most of the order of 10 6 yr [119]. Finally, in cool WDs, when the electron degenerate CO core is in the liquid phase, the slow gravitational settling of 22 Ne lenghtens the cooling times of faint WDs from high metallicity progenitors by up to approximately 1 Gyr, about 10% or more of the total cooling time of the models [123][124][125][126][127]. This element is the most abundant impurity present in the CO core, and originates from 14 N during the core He-burning phase. Owing to the two additional neutrons hosted by the 22 Ne nucleus relative to 12 C and 16 O ions, this element experiences a net downward gravitational force and slow settling in the liquid region of the core (inhibited in the solid regions, due to the expected sudden increase of viscosity). This regime is much denser than the case treated by the Burgers equations, and the 22 Ne diffusion velocity w22 Ne is determined from a diffusion coefficient D22 Ne estimated from molecular dynamics calculations for a strongly coupled plasma [128]: where m p is the proton mass and g is the local acceleration of gravity.   Here, Z 5/3 is an average over the ion charges, T is the temperature and the electron sphere radius a e is equal to (3/4π n e ) 1/3 with n e =Zn the electron density and n the ion density.
In addition, with ω p the plasma frequency ω p = 4π e 2Z2 n M 1/2 , withM denoting the average mass of the ions. The diffusion of 22 Ne increases the energy budget of the WD (hence increases the cooling times) through the contribution (∂U/∂μ) T,V (∂μ/∂t) to the gravitational energy generation coefficient g in the equations of stellar structure (see equation (2.6)).
Moving to stars with initial masses in the range approximately 1.5-3.0M , which display vanishing convective surface layers during the MS, calculations with diffusion show the development of surface abundance variations already during the pre-MS [129]. During the MS phase, they display large surface abundance variations (like He and Ca underabundances, iron-peak element overabundances); moreover, local enhancements of iron and nickel abundances below the thin surface convective layers lead to extra convective zones that, in some cases, can trigger pulsations through the iron-induced κ-mechanism [130].

Inhibition of the efficiency of atomic diffusion
As already mentioned, diffusion is an element transport mechanism that comes from basic physics principles, and should be efficient in stars. Helioseismic observations tell us that diffusion is efficient in the Sun; its inclusion in stellar models improves the match of the inferred sound speed profile, as well as the present depth of the convective envelope and its He abundance [131,132]. The very efficient diffusion of metals from the WD envelopes is also confirmed by observations, whereby for the vast majority of WDs the chemical composition of the surface is either pure H or pure He. 16 Also the diffusion of Ne in the CO liquid cores of WDs is indirectly favoured by the comparison of TO and WD ages for the old metal-rich open cluster NGC6791. The effect of Ne diffusion increases the age obtained from the WD luminosity function, ensuring agreement with the TO age [126].
However, existing spectroscopic measurements of surface element abundances (e.g. Fe) in a sample of GGCs contradict the predictions of stellar models that include uninhibited diffusion [120,[133][134][135][136]. In a star cluster, the initial mass of the stars evolving in post-MS stages is essentially constant, and measurements of chemical abundances at the TO and the base of the RGB should display a difference, due to the effect of diffusion on evolutionary tracks discussed before, which is maximized at the TO and essentially disappears along the RGB. For example, in the case of NGC6397, with [Fe/H]∼ −2.1 along the RGB, the TO Fe abundance should be approximately 0.3 dex lower, and the Ca abundance approximately 0.1 dex larger than on the RGB according to models calculated with diffusion; observations show instead TO Fe and Ca abundances approximately 0.15 and approximately 0.1 dex lower than on the RGB, respectively [137].
This points to a reduction of the effect of diffusion on the surface abundances by some competing mechanism. Nothing of course can be said about the efficiency of diffusion in the inner layers. 17 The same qualitative result is found in more metal-rich clusters like NGC6752 and M4, and also on the HB of a number of clusters, as shown by figure 19. In this figure, one can clearly see that for HB stars hosted by three metal-poor GGCs, the surface Fe abundance increases abruptly when T eff goes above approximately 11 000 K, while theory predicts Fe enhancements (compared with the initial value [Fe/H] = −2.3 dex typical of these clusters) at much lower temperatures.
Also comparisons of models in the mass range between approximately 1.5 and 3.0M , calculated including diffusion with surface abundances measured in A and F stars, show that diffusion must be moderated by some competing mechanism [129,139,140]. For example, in samples of A and F stars belonging to the Hyades, Pleiades and Coma Berenices Galactic clusters, observed abundances of elements like Na, Fe and Ni display a general pattern and star-to-star scatters that point to diffusion being moderated to different degrees by some competing process [141].
Finally, the observed constant surface Li abundance observed in field halo stars with T eff above approximately 5500-6000 K and [Fe/H] below approximately −1.5 dex-the so-called Spite plateau [142]-is also problematic for models including diffusion, because they predict for these stars a progressive decrease of surface Li with increasing T eff , which is not observed [143,144].

Rotational mixing
As we will see in §8, rotational mixing can, in principle, counteract the effect of diffusion, and this has been invoked to explain the delayed (in terms of T eff ) onset of diffusion in hot HB stars [145]. This inference is related to the distribution of observed projected rotational velocities for the hot HB stars of figure 19, as shown in the lower panel of the same figure. All observed stars with T eff above approximately 11 000 K that show spectroscopically the signature of diffusion are very slow rotators, whereas at lower T eff rotation rates are on average higher. Estimates (not based on fully consistent stellar evolution models) of the competing effects of meridional circulation and diffusion, show that moderation/inhibition of diffusion below T eff ∼ 11 000 K by rotational mixing is possible [145]. Interestingly, in support of this idea figure 19 shows three stars with T eff > 11 000 K and large projected rotation velocities with no enhancement of surface Fe, at odds with the other slow-rotating objects at the same temperatures. The link between (partial) inhibition of diffusion and rotation is however probably not justified for A and F stars, which often show slow rotation velocities and surface abundances clearly affected by a reduced efficiency of diffusion, compared with the model predictions. 16 Metals observed in a minority of WD spectra are ascribed to the effect of accretion. 17 If diffusion is inhibited only from the outer convective layers of low-mass stars, the effect on the interiors can still affect the age determination of clusters and field stars [138].

Mass loss
Another process that can moderate diffusion in the stellar outer layers is mass loss [122,146,147]. Assuming that mass loss (hereafter ML) is spherical and unseparated (chemical composition of the wind equal to the photospheric composition), its effect is simply to peel off the outer layers of a model. Simple mass conservation constraints, coupled to the fact that the structure of the star is unaffected by the amount of mass lost between two consecutive computational time steps, translate into the appearance of an outward flowing interior velocity due to the wind [148] given by v w (r) = −Ṁ 4π r 2 ρ m r M * , (6.20) where ρ is the local density at radial distance r from the centre, m r is the mass interior to r, M * is the total mass andṀ (negative) is the mass loss rate. Figure 18 displays the effect of ML on the surface abundances of a typical model at the TO of a metalpoor GGC. For elements affected mainly by gravitational settling, the effect of ML is to increase their surface abundance compared with the pure diffusion case. The reason is that with the chosen mass loss rates, below the convective envelope v w counteracts the settling velocity, hence the degree of depletion is reduced. Larger negative values ofṀ imply larger positive v w and less diffusion from the envelope. In the case of elements supported by radiative levitation in some layers below the convective envelope, the effect of ML on the surface abundances is more complex, and this is due to the interplay between v w and the depth at which g rad is larger than the local gravitational acceleration g. For example, let us assume an ML rateṀ = −10 −13 and an age of 10 Gyr; if v w is larger than the settling velocity down to M = 10 −3 M * from the surface, and in those layers g rad,i > g for a generic element i, after 10 Gyr the wind will have advected to the surface matter at that location, and the surface abundance of element i will be enhanced compared with the initial value.
This complex behaviour of the elements supported by radiative levitation can be seen in figure 18. ML rates of the order of ≈ 10 −12 M yr −1 are required to explain the Spite plateau and the behaviour of surface abundances in the metal-poor GGC NGC6397 [122,146]. These rates are unlikely, being much larger than the solar current mass loss rateṀ = −2 × 10 −14 M yr −1 .
Rates of the order of ≈ 10 −13 M yr −1 are invoked to explain at least some of the surface abundance patterns seen in A and F stars [148].

Thermohaline mixing
The modifications of the chemical abundance profiles caused by diffusion lead to variations of the local mean molecular weight. For example, He always sinks, introducing a stabilizing contribution to the local μ-gradient. On the other hand, heavy elements supported by radiative levitation will lead to a destabilizing μ-gradient (increasing outwards). If this second effect prevails, these layers are subject to thermohaline mixing (see §5), which will tend to erase the abundance changes, and hence to moderate the effect of diffusion.
Published calculations of models for A and F stars, including both diffusion and thermohaline mixing, show that thermohaline mixing (employing the diffusion coefficient derived by Brown et al. [149]) plays an important role in moderating the effect of diffusion on the surface abundances, although the surface abundances of the models do not yet match observations [140]. Both processes have been included also in calculations to explain the observed population of carbon-enhanced metal-poor stellar stars [150].

Generic turbulence
A pragmatic solution adopted in several investigations is to add by hand some turbulence able to counteract the efficiency of atomic diffusion in the outer layers. Turbulence here means simply an additional ad hoc term in the chemical evolution equations that acts towards suppressing the development of chemical abundance gradients, not explicitly connected to a specific physical mechanism. Equation abundance X i of element i, is therefore modified to include an additional term (6.21) where D turb is the chosen turbulent diffusion coefficient.
The literature presents at least three different choices for D turb . To date, the only formulation employed to interpret several different sets of data is the following [118,139,151]: is the He diffusion coefficient 18 at a certain reference depth in the approximation of trace amount in an ionized hydrogen plasma. The subscript 0 indicates the reference depth in terms of either a reference density ρ 0 or, more often, a reference temperature T 0 . In this latter case, ρ 0 is the density at the layer where T = T 0 . The turbulent diffusion coefficient D turb is f times the He diffusion coefficient at the reference depth 0, and varying as ρ n .
A choice of parameters f = 400, log(T 0 ) = 5.5, n = −3 (denoted as T5.5D400-3) allows to match the abundances of A and F stars, while a T6.0D400-3 or T6.09D400-3 choice of parameters in equation (6.22) reproduces the flatness of the Li Spite plateau [118]. Spectroscopy of stars in NGC6397, NGC6752 and M4, three GGCs spanning a range of about 1 dex in initial metallicity, suggests choices between T6.0D400-3 and T6.2D400-3 to match the observations [133,135,152]. Interestingly, for the approximately 4 Gyr old, solar metallicity open cluster M67, predictions from stellar models with uninhibited diffusion are generally consistent with observed abundances [153], although the predicted abundance variations are small.
Two additional forms for D turb have been proposed, but not widely applied, to comparisons with spectroscopic observations. The first one applied to low-mass star models with convective envelopes is (6.23) with f = 15 and n = 1.5 [154], where BCZ denote quantities at the lower boundary of the convective envelope and M tot is the total mass of the model. Owing to the steep power-law dependence of the density ratio, D turb becomes negligible in the nuclear-burning regions, ensuring that the assumed turbulence will not affect the inner chemical profiles resulting from nucleosynthesis and gravitational settling. The other proposed parametrization is proportional to the radiative viscosity [155] D turb,visc = f 4aT 4 15cκρ 2 , (6.24) where f is a free parameter that is constrained in the range f = 1 +2.0 −0.2 by helioseismic observations, and spectroscopy of Hyades stars and OB stars (in this latter case assuming no competing effect from ML). Figure 20 compares the different formulations for D turb discussed above, for the outer layers of a 0.8M , [Fe/H] = −1.3 model, in the latter phase of its MS evolution. In the case of D turb,Rich , we display both the T6.0D400-3 and T6.2D400-3 results, the ones that allow a match to results from spectroscopic observations, as discussed before. Notice how when moving the reference temperature from log(T 0 ) = 6.0 to log(T 0 ) = 6.2, the entire profile of D turb,Rich shifts to deeper layers, thus counteracting the efficiency of diffusion in deeper stellar regions. Also, D turb , VdB (displayed with the choice of parameters recommended in [154]) appears to be very similar to D turb,Rich for the T6.2D400-3 choice. On the other hand, both absolute values and trends of D turb,visc (for f = 1) with mass depth are very different from the other cases.
For the sake of comparison, we also show the diffusion coefficient of He (that is never supported by radiative levitation) at the edge of the convective envelope and below. All D turb choices displayed in the 18 The He settling velocity w He can be related to a diffusion coefficient D(He) through where m He is the mass of the He nucleus. figure are able to counteract the diffusion of He from the bottom of the convective zone. Even D turb,visc , that is the smallest turbulent diffusion coefficient in the figure, matches the one for He, hence is able to inhibit its settling below the convective envelope.

Phase separation in white dwarfs
Besides convection in the non-degenerate envelopes and the diffusion of 22 Ne in the liquid core, another important process that redistributes the chemical abundances in cold WDs is the so-called 'chemical separation' upon crystallization [156]. In a nutshell, when a layer with given C and O abundances in the liquid phase crystallizes (we assume for the moment the core being made of just C and O, given that the sum C+O is always more than approximately 95% by mass in any realistic WD model), the equilibrium composition in the solid phase-given by the phase diagram of the CO mixture-is very likely to be different. This will cause a gradual change of the overall chemical profiles within the CO core that impacts the energy budget and cooling times, as detailed below. Figure 21b displays a phase diagram for the CO mixture [158], with the regions of liquid and solid phase labelled. For heuristic purposes, let us assume a C mass fraction X C = 0.50 uniformly throughout the core (hence the oxygen mass fraction X O is also approx. 0.50). When the mixture starts to crystallize at the centre (the value of the Coulomb parameter Γ = 180 for the onset of crystallization is reached first in the centre, because of higher densities), the chemical composition in the solid phase is determined as follows. One needs to draw a vertical line in the phase diagram of figure 21 with the horizontal coordinate equal to 0.50, which runs through the region belonging to the liquid phase until it intersects the upper line describing the phase diagram. The vertical coordinate of the intersection point gives the crystallization temperature of the WD centre (corresponding to T ∼ 1.3T C , where T C denotes the crystallization temperature of a pure C mixture). From this point, one has to draw a horizontal line that will intersect the lower segment of the phase diagram in correspondence with a carbon abundance X C ∼ 0.30. This is the equilibrium abundance of carbon in the solid phase (the corresponding oxygen abundance is X O = 1 − X C ).
Given that X C in the now-crystallized centre is lower than the initial value, conservation of mass requires that the carbon abundance in the liquid phase at the crystallization boundary is increased (hence X O is decreased) with respect to the original value. This means that right above the crystallized boundary the molecular weight is now lower than in the overlying layers still in the liquid phase (where the ratio X O /X C is higher). An increase in molecular weight with increase in distance from the centre in the liquid phase causes an instability to develop, and the resulting fast-mixed region extends outwards in mass   [157]. (b) Phase diagram of the CO mixture adopted to calculate the profiles in the left panel [158]. T C denotes the crystallization temperature of a pure C composition, X C the C mass fraction.
as long as the new uniform average X C value is higher than the abundance in the next, unperturbed layer (in this case, it will reach the edge of the CO core) [159,160]. After mixing is completed, the liquid phase will have an overall enhanced X C abundance compared to the value before the crystallization of the central layer.
Let us now suppose that the new value of X C at the boundary of the solid core is equal to 0.55 when this layer crystallizes-at a lower temperature than the core because of lower density. The abundance in the solid phase can be derived in the same way as before, and it is now equal to X C ∼ 0.35. This implies again that right outside the newly crystallized layer the abundance X C must be higher (and X O lower) than in the overlying layer. The instability in the liquid phase ensues again and the cycle is repeated (mixing in the liquid phase eventually stopping when the carbon abundance of the newly crystallized layer becomes lower than or equal to the overlying layers still in the liquid phase) until the whole degenerate core is crystallized. The final profile of X C and X O after crystallization is completed is no longer homogeneous; X O will display central values higher than in the liquid phase which decrease from the centre outwards, while the opposite is true for X C .
This variation of the X O /X C profile (hence the local value of the molecular weight μ) in the WD core during crystallization due to the CO phase diagram has an important impact on the WD cooling times, because of the term (∂U/∂μ) T,v (∂μ/∂t) in g [157]. This means that more energy is available to be released and cooling times are longer. Depending on the thickness and chemical composition of the non-degenerate layers surrounding the CO core (that regulate the energy release), delays induced by chemical separation upon crystallization can reach approximately 10% for the coolest objects. Figure 21a displays the evolution of the O-profile for a realistic 0.61M WD model [157]. The first stage of this temporal sequence is represented by the profile in the electron degenerate core at the beginning of the thermal pulse phase of the AGB progenitor, built during the core He-burning. One can notice the local maximum in X O is about 0.25M away from the centre. The associated local maximum of the molecular weight μ causes an instability that homogenizes the internal layers during the liquid phase of the cooling [157], generating the profile that is maintained until the start of core crystallization. The final very different X O profile is attained when the whole core is crystallized, and it is produced by the chemical separation upon crystallization.
The most recent re-evaluation of the phase diagram for a binary CO mixture is displayed in figure 22 [161]. Compared to the widely used phase diagram of figure 21, it causes slightly smaller abundance variations for a given initial chemical profile in the liquid phase.
Within the CO core there are also small amounts of other metals-so-called minor species, with mass fractions of at most the order of the initial progenitor metallicity-that have been either processed through the previous burning phases (i.e. Ne), or are unchanged since the formation of the progenitor  Figure 22. The most recent re-evaluation of the phase diagram for a CO mixture [161].
(i.e. Fe). Owing to the extreme complexity of calculating a multi-component phase diagram, a ternary CONe mixture is often assumed to behave as an effective binary mixture composed of neon (or iron) plus an element of average charge Z = 7 determined from the C and O abundances. Ne is important, given that 22 Ne is the most abundant species after C and O, with a mass fraction X22 Ne ∼ Z, where Z is the initial metallicity of the progenitor (Z ∼ 0.02 for the solar metallicity). A detailed phase diagram for a three-component mixture CONe shows that the effect of Ne on the final WD cooling times is negligible, when the separation of carbon and oxygen is accounted for [162].

Rotation and rotational element transport mechanisms
Basic considerations about the angular momentum evolution in a contracting protostellar cloud, observations of the solar magnetic field and stellar spectroscopy (and now also asteroseismology) dictate/show that stars do rotate. Observations of MS O-and B-type stars in the Galaxy and the Magellanic Clouds reveal average projected rotational velocities of the order of 150 km s −1 [163][164][165][166], with values up to 350-400 km s −1 . Average rotational velocities of the order of 150 km s −1 are observed also in MS A-and F-type stars [167], fast decreasing when moving to later spectral types [168], while giant stars display typically slow projected rotational velocities below approximately 10 km s −1 [169]. All fast rotators among O-stars show surface He-enrichments not predicted by standard non-rotating models [170]. Also, high rotational velocities along the MS are associated with enhancements of surface N that cannot be reproduced by non-rotating stellar models [171,172].
For a long time, rotation has not been an ingredient of standard stellar models, due to the increase in complexity and uncertainty-free parameters-related to the inclusion of rotation. Moreover, the basic principle that the explanation relying on the smallest number of hypotheses is the one to be preferred, coupled to the many successes of non-rotating stellar models, means that rotation has been generally considered only a second-order effect. It is, however, clear that it can have important effects on the structure and evolution of stars, mainly through its effect on the evolution of the chemical stratification.

One-dimensional modelling of rotating stars
The physical basis for the inclusion of rotation in detailed stellar evolution computations was laid down 50 years ago [173,174] and more recently hydrodynamical treatments of the problem have appeared [175][176][177]. Detailed stellar evolution calculations and comprehensive libraries of stellar models are, however, still (and for the foreseeable future) possible only with standard 1D modelling, and the basic (i) Roche approximation, i.e. the gravitational potential Ψ is the same as if the total mass of the star were concentrated at the centre. (ii) The angular velocity Ω and chemical composition are constant along isobars, e.g. surfaces of constant pressure (shellular rotation). This follows the assumption that turbulence is anisotropic, with a stronger transport in the horizontal (tangential to an isobar) direction than in the vertical (perpendicular to an isobar) one [178]. The angular velocity in this case depends very weakly on the colatitude, hence Ω ≡ Ω(r), meaning that it varies according only to the radial coordinate of the isobars. (iii) The mean radius r P of an isobar is defined by V P = (4π/3)r 3 P , where V P is the volume inside the isobar.
With these assumptions, the set of equations of stellar structure for a rotating star can be written in 1D as [179,180] ∂r with ∇ P being the appropriate temperature gradient d ln(T)/d ln(P). In the case of radiative transport, When comparing these equations to equations (2.1)-(2.4), one can notice that they are formally identical, apart from the form factors f P and f T . The only difference is the interpretation of the variables. In the case of rotating models, r P is the mean radius of an isobar defined by the value P of the pressure,ρ andT are the volume-averaged density and temperature between two contiguous isobars (the difference with averages on the isobars is negligible if the mass grid of the models is dense enough), and L P and m P are the mass and luminosity inside a given isobar. Calculations of the energy generation coefficients, adiabatic gradient and opacity make use ofρ andT. The equation of state is expressed in terms of P = P(ρ,T) and the chemical composition on the isobar P.
Analogous to the case of non-rotating models the system of equations is solved considering the 'Lagrangian' independent variable m P , with r P , L P , P,T as unknowns.
The form factors f P and f T are defined, respectively, as (8.4) The quantities g and g −1 are average values of the gravity g over an isobar P with surface area S P . For a generic variable f this average is defined as where dσ is an infinitesimal element of the isobar surface P. At the non-rotation limit, f P and f T converge to unity and the equations are reduced to their non-rotating form. To solve equations (8.1), one needs to evaluate f P and f T , and this requires the calculation of the surfaces of isobars, e.g. surfaces where Ψ P is constant, with Ψ P = −Φ + 1 2 Ω 2 r 2 sin 2 θ, (8.6) where Φ is the Roche gravitational potential, r is the radial distance from the centre and θ is the colatitude (θ = 0 at the poles; see e.g. [181] for an example of how to calculate integrals (8.5) for a given Ψ P ). It is through the calculation of f P and f T that the angular velocity profile enters the equation of stellar evolution. It is worth noticing that all this formalism works also in the case of a 'conservative' potential (for example, the simple case of solid body rotation), e.g. if the centrifugal acceleration can be derived from a potential [182,183], which is not the case for shellular rotation. A potential of this type is used, for example, in [182] In this case, instead of isobaric surfaces one reads equipotential surfaces, and the equations (including f P and f T ) are the same as the shellular case.

Chemical element and angular momentum transport
A tricky issue for modelling rotating stars with 1D stellar models is how to describe the transport of angular momentum-that determines the evolution of Ω-and the chemical mixing associated with rotation. These are two facets of the same problem, for rotation triggers hydrodynamical instabilities and large-scale motions of the gas in radiative regions, which result in transport of both angular momentum and chemical elements. In-depth discussions about rotation-driven instabilities can be found in [180,184,185]. Here, we just give a brief overview, focussing on the actual implementation in stellar evolution calculations. The Eddington-Sweet meridional circulation is one of the major instabilities caused by rotation. In simple terms, we can compare a non-rotating star in radiative equilibrium with its solid body rotating counterpart. The equipotential surfaces of the non-rotating star are spherical, while in the case of a solid body rotator they are rotational ellipsoids, and two contiguous equipotential surfaces will diverge in distance from each other at the equator. Given that the effective gravity g is proportional to the gradient of the potential (that is normal to the equipotential surfaces), it will vary with latitude on an equipotential surface. As a consequence, the temperature will be hotter at the poles and cooler at the equator-as demonstrated by von Zeipel almost a century ago, the energy flux is proportional to the local value of g (the von Zeipel theorem)-preventing the star from maintaining hydrostatic equilibrium.
The solution to this paradox is to invoke large-scale mass motions that transport energy, the socalled meridional circulation, moving material inwards from the equator and upwards along the rotational axis towards the poles. The timescale for this mixing process-the 'Eddington-Sweet' timescale-was estimated to be where L, M and R are the stellar luminosity, mass and radius, respectively, Ω the angular velocity and the first term is the Kelvin-Helmholtz thermal timescale. These early estimates of t ES were much shorter than the MS lifetime of stars, even for modest rotation rates, hence rotating stars should be fully mixed, contradicting observational data. The presence of chemical abundance stratifications (μ-gradients) in the interior of rotating stars can, however, increase t ES considerably compared to the early estimates. A widely used model of the meridional circulation included in stellar evolution codes (that will be used later when discussing the implementation of rotational mixing) develops the circulation velocity vector U into two components [178]: U = U 2 (r)P 2 (cos(θ ))e r + V 2 (r) dP 2 (θ) dθ e θ , (8.9) where r is the radial coordinate, θ is the colatitude and P 2 is the Legendre polynomial of order 2. The radial components of these velocities, which will be used to treat chemical element transport, are related through 1 r loss) will generate-starting, for example, from solid body rotation, usually assumed for pre-MS fully convective stars-a variation of angular rotation velocity with depth in radiative regions. Hence a 'shear' develops between neighbouring layers that leads to instabilities. This stems from the fact that the minimum energy state of a rotating fluid is solid body rotation, and if the star develops differential rotation, it is possible to extract energy by homogenizing the velocities through transport of material.
The strong thermal and molecular weight radial stratification in radiative zones tends to oppose the homogenization of the rotational velocities, while it is reasonable to assume that no restoring forces oppose horizontal displacements. As a consequence, horizontal shear (along an isobar) is expected to generate a strong turbulence on short dynamical timescales, justifying the assumption of shellular rotation. Eventually, thermal diffusion and horizontal shear can reduce the stabilizing effect of the vertical (radial) thermal and chemical stratification and induce element and angular momentum transport if the Richardson number Ri satisfies this criterion: where u is the velocity of the fluid elements and z designates the vertical direction. 20 The Solberg-Høiland instability, plus additional instabilities (i.e. Goldreich-Schubert-Fricke, ABCD instabilities) that can develop when equipotentials and isobars do not coincide-so-called baroclinic instabilities [180]-do potentially contribute to the transport of angular momentum and chemicals. 21 Some of these instabilities are implemented in calculations of rotating stellar models (see §8.2.2), even though their treatment is considered to be much more uncertain than the case of meridional circulation and shear instabilities [184,185].
In the case of radiatively driven winds, the von Zeipel result has important consequences on the mass loss rates to be employed in the calculation of rotating stellar models, because in general it causes an increase in the mass loss efficiency compared to the case of a non-rotating counterpart. Various prescriptions of the mass loss enhancement due to rotation can be found in the literature. As an example, we report the prescription employed in the MESA code [58]: (8.12) with ζ =0.43, and Ω 2 crit = (1 − L/L Edd )GM/R 3 and L Edd = 4π cGM/κ averaged over a certain optical depth range.

Advective/diffusive implementation
Meridional circulation and shear instability are considered as the main mechanisms for angular momentum transport and rotational mixing in a number of stellar evolution codes, e.g. STAREVOL [187,188], CESTAM [189], the Geneva stellar evolution codes [179,190] and the FRANEC code [191], that implement advective+diffusive transport of angular momentum and diffusive rotational chemical mixing.
The transport equation for a chemical element with mass fraction X i is written as 22 where D chem is the sum of the vertical shear diffusion coefficient D shear and the effective diffusion coefficient D eff , which accounts for the combined effect of the strong horizontal shear diffusion D h and of the meridional currents: 14) The equation for the transport of angular momentum is written as [192] ∂ ∂t (r 2 Ω) where D ang denotes the diffusion coefficient for angular momentum. The second term on the right-hand side of equation (8.15) is a diffusion term, similar in its form to equation (8.13), while the first term is an advective term, modelling the transport by a velocity current. Equation (8.13) does not contain an advective term, following the author of [193] who showed that the combined effect of turbulence and circulation currents is equivalent to a diffusion process for the element transport.
There are various expressions in the literature for U 2 . The complete description by Maeder & Zahn [194] in the case of shellular rotation provides The barred symbols indicate averages over the isobar; E Ω and E μ are complicated terms that depend on the angular velocity profile and mean molecular weight fluctuations on an isobar (see [180] for the complete formulae and the derivation); Θ describes the density fluctuation on an isobar; and M here denotes the reduced mass where ρ m is the mean density inside the considered isobar surface, and L is the luminosity within the same surface.
In the case of convective regions, the transport of angular momentum is very uncertain, as the interaction between convection and meridional currents is not well understood. Traditionally, one can choose between the following two limiting cases. If convection inhibits meridional currents, angular momentum is redistributed very efficiently by convection, as in the case of chemical elements, and the result will be solid body rotation, as in the solar convective envelope. If meridional currents dominateone can hypothesize that this is the case applicable to large, rarefied RGB envelopes, where convective elements may collide elastically rather than inelastically as in the solid body case-it is the specific angular momentum that is expected to be uniform. Hydrodynamical simulations have shown that the extended deep convective envelopes of red giant stars are likely to undergo radial differential rotation, with an angular velocity profile of the form Ω(r) ∝ r −0.5 [195]. In general, solid body rotation is usually prescribed in convective regions. Angular momentum losses due to stellar winds need also to be accounted for. In the assumption that the mass loss is spherically symmetric, the rate of angular momentum loss during a given time step will be approximately equal to the average specific angular momentum at the surface of the star multiplied by the assumed mass loss rate.
To discuss the effect of rotation and rotational mixing on stellar evolution models, we consider as an example models from [196], calculated with the Geneva code. Figure 23 compares the HRD of 9, 15 and 32M solar initial chemical composition models with and without rotation, from the zero age main sequence (ZAMS) to the end of core C-burning. In these calculations, D h has been taken as from [178], assuming c h = 1 and α = 1 2 (d ln(r 2 Ω)/d ln r). We express D shear as from [197], where f energ = 1 (the fraction of the excess energy in the shear that contributes to mixing). The models are initialized as solid bodies at the ZAMS, and then evolved according to the prescriptions for the angular momentum and chemical transport described above. The prescribed initial equatorial rotation velocity is 0.4 times the critical velocity (when the centrifugal acceleration in the equatorial plane exactly compensates the gravitational acceleration) of the corresponding stellar mass. This is equal to values between approximately 250 and 300 km s −1 for the three displayed masses.
The effect of rotation on the tracks (that of course depends on initial velocities) is striking. Models calculated with rotation evolve to higher luminosities during the MS and stay more luminous also during the following evolutionary phases. This is mainly due to the consequences of the chemical element transport associated with rotation. In general, rotational mixing works in the direction of erasing chemical gradients within the star, and hence during the MS, this causes a continuous slow ingestion of fresh hydrogen into the H-depleted convective core, as well as a slow transport towards the surface of elements whose abundances are increased by H-burning (e.g. He and 14 N). The effect of the ingestion of  (c) ( d) Figure 24. Internal profiles of H and N abundances, angular velocity and chemical diffusion coefficients for shear and meridional circulation, respectively, as a function of mass, within 15M solar initial composition models from [191], with the labelled initial rotation rates. The different panels refer to two MS stages, when the central H mass fraction is equal to 50% (a,c) and 10% (b,d), respectively (courtesy of M. Limongi).
H in the convective core is to increase the MS lifetime compared to the non-rotating counterpart, and to produce slightly larger He-cores at the end of core H-burning.
To this purpose, figure 24 displays the internal profiles (taken during the MS when central H is reduced to-from (a,c) to (b,d)-50% and 10% by mass, respectively) of H and N abundances, angular velocity and chemical diffusion coefficients for shear and meridional circulation in solar initial composition, 15M rotating and non-rotating evolutionary models from [191], similar to the Geneva calculations. One can see very clearly the effect of rotational mixing, which smooths out gradients in the N-and H-abundance profiles. Also, shear mixing dominates in the external layers, where the gradient of angular velocity gets progressively steeper, while meridional circulation is the more efficient chemical  transport process close to the edge of the convective core. Notice also the flat angular velocity profile in the convective core, due to the solid body assumption.
The different paths in the HRD followed by rotating and non-rotating models together with the different lifetimes lead also to different mass loss histories and hence to different total masses at the end of the MS and during the post-MS phases. This impacts, for example (for a fixed convective mixing scheme), the occurrence of 'loops' in the HRD, as the blue or red location of a model during the giant phase depends on thickness of the H-rich envelope (thick enough envelopes keep the models at the red side of the HR, while when their masses decrease below a threshold value, the models move to the blue). During the He-burning phase, the slow ingestion of fresh He in the core lowers the final C/O ratio and increases the mass of the CO core [191], affecting the pre-supernova structure of the massive models. Figure 25a shows the evolution of the angular velocity Ω taken at the surface and at the centre of the 15M models. After the ZAMS model with enforced solid body rotation is left to evolve, there is a readjustment of the rotational profile at the very beginning of the evolution, until the equilibrium rotational profile is reached, and the angular velocity starts to evolve under the action of the transport mechanisms. The general trend is to transfer angular momentum from the contracting core towards the external layers, but this is counterbalanced by the mass loss that removes angular momentum from the surface and eventual expansions of the convective envelopes that slow down the surface.
For the model shown in figure 25, the net effect is a constant decrease of the surface angular velocity along the whole evolution, whereas the centre of the star displays a moderate increase of angular velocity during the MS, followed by a plateau and a sharp increase during the giant phase. The models are never very far from spherical symmetry, as shown in figure 25b, which displays the ratio of the polar to the equatorial radius along the whole evolution.
The effect of element transport on the surface abundances (in mass fractions) of some key elements (He, C, N, O) for the same 15M evolution is displayed in figure 26. Notice that in the non-rotating models the abundances change only due to the dredge-up during the red giant phase. Rotating models display instead abundance changes already during the MS, due to rotational mixing that tends to erase chemical gradients. This explains the increase of 14 N and He, whose abundance increase in the central regions due to CNO H-burning, and the depletion of 12 C and 16 N, caused by the decreased O and N equilibrium abundances, compared to the initial scaled solar values. The abundances remain almost constant during the fast transition to the giant phase, and then change again due to the dredge-up.
The detailed behaviour of the evolution of the angular momentum and chemical abundance profiles is, however, strongly dependent on the precise choice of the diffusion coefficients in equations (8.13) and (8.15). We can see, for example, that the coefficient employed in these calculations contain parameters like c h and f energ that are set to fixed constant values not derived from first principles.  It is extremely interesting to analyse also the case of the 1M rotating (ZAMS equatorial velocity equal to 50 km s −1 ) and non-rotating models for the same initial chemical composition. Both calculations account also for atomic diffusion, but without including radiative levitation, whose effect is practically negligible at this metallicity. Figure 27 compares their HRDs, 23 which are virtually identical along the MS, while the rotating SGB is slightly brighter and the following RGB phase slightly hotter. The MS lifetime of the rotating model is just approximately 4% longer than the non-rotating case. then increase when convection deepens, before the signature of the first dredge-up can be seen for C, N and He. When rotation is included, the abundances stay constant along the MS, almost equal to the initial values, showing that rotational mixing (for the chosen initial rotational velocity at this metallicity) strongly inhibits the effect of diffusion from the convective envelopes. Along the RGB of the rotating models displays stronger enhancements of N and He, and a larger depletion of C compared to the nonrotating counterparts. There are also alternative expressions for D shear and D h , in addition to the ones employed in the calculations we are discussing, namely: D shear from [198] D shear = f energ H P gδ with K, f energ and ϕ as in [197]. D h from [199] D h = Ar(rΩV 2 (r)|2V 2 (r) − αU 2 (r)|) 1/3 , (8.21) with α as in [178] and A = 0.002. D h from [200] D h = β 10 with α as in [178] and β = 1.5 × 10 −6 .
The effect of using these various combinations of D shear and D h in the Geneva code has been investigated recently [201]. 24 They found that MS lifetimes, the evolution of surface velocities and the angular momentum of the core have a weak dependence on the choice of these diffusion coefficients. The shape of the evolutionary tracks, the surface enrichment (for a fixed initial rotation velocity), the blue-to-red evolution in the HRD and the extension of the blue loops are, however, significantly affected by the choice of D shear and D h .

Diffusive implementation
There is an alternative approach to include transport of chemicals and angular momentum, used in codes like the Yale evolutionary code [203] where ν and D are, respectively, the total turbulent viscosity and the total diffusion coefficient defined as a sum of all diffusion coefficients associated with all the transport processes taken into account. Each of these diffusion coefficients is built as the product of the velocity v and the path length l of the redistribution currents with By |∂r/∂ ln v| is denoted the velocity scale height. The timescale associated with the redistribution over the path scale is simply l 2 /D = l/v. This description is usually employed to include within the same formalism also convection and semiconvection, treated as diffusive processes. Appropriate diffusion coefficients for convective and rotational transports are adopted, often with free parameters to be calibrated against some sets of observations, given that these coefficients arise from order-of-magnitude considerations. Rotational mixing processes include at least meridional circulation, and shear, plus eventually additional rotational instabilities [58,76] not usually included in the codes that employ the advective/diffusive implementation.
The viscosity ν is prescribed by (8.26) and the diffusion coefficient is usually written as the following sum [76]: where f c is one of the efficiency parameters entering the diffusive formalism and it is calibrated on observations like the solar lithium abundance (f c = 0.046 in the models [203]), or the main trend of the observed nitrogen surface abundances with the projected rotational velocity for the nitrogen-enriched fast rotators in the LMC (f c = 0.0228 in the models [204]). An example of these implementations is given in the following, for the meridional circulation [76,203]. The characteristic velocity of these currents in a radiative region is assumed to be , (8.28) with denoting the energy generation rate per unit mass. The counteracting effect of μ-gradients is accounted for according to ψ = (d ln(ρ)/d ln(μ)) P,T , τ KH is the local thermal timescale and f μ is a free parameter that allows to vary the sensitivity of meridional circulation to μ gradients.
The effective value of the meridional circulation velocity used to determine the associated diffusion coefficient is finally calculated as v ES = |v ES 0 | − |v μ | if v ES 0 > v μ , or v ES = 0 otherwise. Figures 29 and 30 compare results for the initial solar composition 15M models discussed in figures 23, 25 and 26, with results of a MESA calculation (that employs this diffusive implementation for the transport of angular momentum and chemicals) for the same mass, solar initial chemical composition and approximately the same initial rotational velocity [205]. The HRD shows that these rotating models are almost equivalent to the non-rotating Geneva calculations. However, the evolution of the surface chemical elements shown in figure 30  quantitative effect is very different from the Geneva results. The enhancement of the surface abundances during the MS is more moderate in the MESA calculations (even accounting for the slightly different initial abundances of some elements), and He is almost unchanged. Even the effect of the dredge-up seems to differ between the calculations.
It is difficult to determine the exact cause of these differences, given that also other elements of the model input physics (apart from the implementation of rotational transports) differ between the two sets of calculations. Very importantly, MESA calculations include also the effect of atomic diffusion (and radiative levitation), however, moderated by some additional turbulence [205]. But taken at face value, this comparison may give an idea of the uncertainties in the outputs of the current generation of rotating models.

Additional effects
There are at least two additional processes that may affect substantially the evolution of rotating stars, for they modify the transport of angular momentum and, directly or indirectly, also the transport of chemical elements. Their implementation in stellar model calculations is still uncertain and they are not generally included in the current generation of stellar models.

Gravity waves
Inside a star, the density changes monotonically with the depth and the molecular weight typically increases in the direction of increasing gravity. When a gas element in this stable stratification is perturbed, the competition between buoyancy and gravity gives rise to an oscillatory motion around the equilibrium position and creates a so-called 'internal gravity wave' (IGW).
These IGWs are expected to be generated by the injection of kinetic energy from a turbulent region into an adjacent stable region, as observed both in two-dimensional (2D) and 3D simulations of convective mixing, e.g. by convective overshooting in a stable region and bulk excitation or excitation by Reynolds stresses inside the convection zone [206][207][208]. The IGWs penetrate into the radiation zone, transporting angular momentum that is deposited where they are dissipated through heat diffusion by photon exchange, which produces an 'attenuation factor' (τ ) proportional to the thermal diffusivity and inversely proportional to the IGW frequency ν and amplitude. It is by shaping the internal rotation profile that IGWs contribute indirectly to the mixing of chemical elements. Angular momentum transport by gravity waves is seldom included in evolutionary stellar model calculations. The current approximate treatment implemented in some stellar evolution calculations [184] expresses the solution of the equations describing the propagation of IGWs in a rotating star in terms of Legendre polynomials. At each point within a radiative region, the total angular momentum 'luminosity' 25 associated with the IGW propagation can be written as where r c denotes the radiation/convection interface, and m represent, respectively, the spherical order and the azimuthal number of the Legendre polynomial, ν the frequency of the IGW when launched from the convective zone, σ = ν − m(Ω(r) − Ω(c)) is the local wave frequency measured in the co-rotating frame with angular velocity Ω(r)-Ω(c) being the angular velocity of the solid body rotating convective zone-and Here, F J (ν, , m) is the mean flux of angular momentum carried by a monochromatic wave with emission frequency ν. The amplitude of the waves is assumed to be [209] F E (ν, ) = ν 2 4π dr ρ 2 r 2 ∂ξ r ∂r with ξ r and [ ( + 1)] 1/2 ξ h being, respectively, the vertical and horizontal displacement wave functions normalized to unit energy flux at the edge of the considered convection zone, v c the convective velocity, Λ = α MLT H P the radial size of a convective element, τ Λ ∼ Λ/v c the convective timescale and h ν = λ min{1, (2ντ Λ ) −3/2 }. The radial (k r ) and horizontal (k h ) wave numbers are related by The local damping rate τ (r, σ , ) can be written as  Figure 31. Evolution with T eff of the central and surface angular velocity, for 1M initial solar chemical composition rotating models, calculated with the Geneva code [196].
where N 2 T is the thermal part of the Brunt-Väisälä frequency and ν v is the vertical shear turbulent viscosity Ri crit (r(dΩ/dr)) 2 35) with N 2 μ denoting the molecular weight stratification part of the Brunt-Väisälä frequency, and ν h the horizontal shear turbulent viscosity, which can be set to D h , and Ri crit = 1 4 . The deposition of angular momentum is then given by the radial derivative of L J . Given that, as a first approximation, only the radial dependency of IGW transport is considered, all quantities required are evaluated from horizontal (on isobars) averages. The angular momentum evolution due only with IGW transport is given by The + or − sign in front of the angular momentum luminosity corresponds to prograde (m > 0) or retrograde (m < 0) waves.
It is clear that in the absence of differential damping for inward-and outward-travelling waves, there will be no angular momentum transport. This is the case, for example, of solid body rotation. However, differential rotation naturally filters the IGWs that propagate within the star. Considering, for example, the realistic case of a low-mass star with interior layers-see e.g. Figure 31 for the 1M solar chemical composition Geneva rotating model discussed before-that rotate faster than the convective envelope. In the inner radiative regions, we have Ω(r) > Ω(c) and for prograde m > 0 waves σ = ν − m(Ω(r) − Ω(c)) is shifted to lower frequencies when travelling from the edge of the convective envelope towards the interior layers. The radial damping length of the waves given approximately by progressively decreases for these prograde waves. The retrograde m < 0 waves are instead boosted to higher frequencies, and their damping length increases, propagating further within the star compared to the prograde waves. As a consequence, prograde waves are absorbed before they can propagate far into more rapidly rotating layers of a star, while retrograde waves pass through, and upon dissipation they deposit their negative angular momentum (they contribute with a negative sign in equation (8.36)), spinning down the rapidly rotating layers. The star then tends to evolve towards solid body rotation. The angular momentum transport associated with IGWs is considered to be a possibility to explain the solar rotation profile, which is much flatter than generally predicted by rotating models that do not include IGWs [210]. Also, recent asteroseismic observations [211,212] show that the core of low-mass RGB stars do not rotate faster than the surface as much as predicted by current rotating stellar models. The angular momentum redistribution by IGWs may provide one way to reconcile theory with observations. Recent multidimensional detailed hydrodynamical simulations of generation and propagation of IGWs within stars, confirm their ability to transport efficiently angular momentum [213,214].

Magnetic fields
Another way to favour the redistribution of angular momentum within a rotating star is to consider the effect of internal magnetic fields. They are generally implemented following the dynamo mechanism presented in [215], which envisages the creation of magnetic fields in the radiative regions of differentially rotating stars at the expense of the shear, due to the so-called Tayler-Spruit instability. 26 The STERN code, which uses a diffusive implementation of element and angular momentum transports, accounts for the effect of magnetic fields as follows [217].
Defining by q = d ln Ω/d ln r the shear, the effective radial viscosity produced by the magnetic field can be written as with η = 7 × 10 11 ln(Λ)T −3/2 denoting the Spitzer magnetic diffusivity and ln(Λ) the Coulomb logarithm (see §6). These viscosities and diffusivities have been included as additional terms to the diffusion coefficients of, respectively, angular momentum and chemicals, in both the STERN [204,217] and KEPLER [218] codes, which use a diffusive implementation of element and angular momentum transports. In this latter code, modifications are applied to the viscosities and diffusivities derived before in the case of semiconvective regions, or where the thermohaline mixing is efficient, which means in regions where the timescale of mixing is much longer than in convective layers. In the case of semiconvection, a dynamo effective viscosity is defined as ν re = ν e0 f (q) and the final expression for the viscosity in the semiconvective region is set to be ν e = √ ν re ν SC , where ν SC = (1/3)H p v conv and v conv is the convective velocity derived from the MLT. The effective diffusion coefficient for the transport of elements is set to (D e + D SC ), with D SC being the diffusion coefficient due to semiconvection. In the thermohaline mixing region D e is set to D e1 , ν e to ν e1 and q min to q 1 . The inclusion of internal magnetic fields decreases the rotation velocity contrast between core and envelope during the MS, approaching near solid body rotation. This affects both the shape of the evolutionary tracks as well as the abundance profiles within the models. Figure 32 displays, as an example, the effect on the surface chemical abundances at selected evolutionary stages in the calculations with the code KEPLER [218], for 15M solar metallicity models with a ZAMS equatorial rotation velocity of 200 km s −1 .
Internal magnetic fields generated by the Tayler-Spruit instability have been also included in the Geneva code [219]-which implements the advective/diffusive formalism for rotational transports-in a slightly different form. Denoting by x the ratio where ω A is the Alfvén frequency of a magnetic field of intensity B, the solution of the equation provides the unknown quantity ω A . The diffusion coefficient for the vertical transport of angular momentum is then given by If this condition is not realized, a formalism for the treatment of low rotation rates-expanding upon work in [215]-is implemented. We notice that, in this treatment, the Brunt-Väisälä frequency is rewritten as A set of 15M solar metallicity models calculated with this implementation [219] also predict an almost solid body rotation during the MS, with obviously a negligible element transport by shear mixing. The transport of elements due to the magnetic fields is also negligible, while the meridional circulation is strongly enhanced by the flat rotational profile. The net effect during the MS is an enhancement of the surface abundance variations compared to rotating models without internal magnetic fields. This is different from the result of figure 32 for a KEPLER model with the same mass, initial chemical composition and very similar initial rotation rate. In this latter case, the variation of the surface abundances during the MS is almost negligible compared to the case without magnetic fields.
Finally, we mention very briefly the effect of magnetic braking due to magnetized winds. 27 Winds with magnetic fields-in the case of surface magnetic fields, whatever their origin is-exert a braking torque that is significantly larger than that for non-magnetic cases [220,221]. The reason is that the magnetic field connects the mass lost from the surface of the star to the envelope, and when the wind finally decouples from the magnetic field, it has the same angular velocity as the surface but a much larger moment of inertia. This increases the amount of angular momentum lost, compared to the case with no magnetic fields, and it is accepted as the reason for the slow rotation rate of low-mass MS stars. The effect on the chemical element transport is indirect, through the variation of the star rotational profile. Surface magnetic fields are observed mainly in low-mass stars with convective envelopes, where they are expected to be generated by a dynamo mechanism, and this type of magnetic braking is implemented mainly in models of low-mass stars. There are various prescriptions in the literature that require the calibration of free parameters on observations, given the lack of knowledge about the surface magnetic fields [222]. For example, the rate of angular momentum loss due to magnetized winds employed in [222] is a variation of [223,224] for Ω < Ω sat , (8.55) where K W = 2 × 10 48 in cgs units to reproduce the solar case and Ω sat is a free parameter.
9. An example of possible synergy among several element transport processes After the description of all major element transport mechanisms included in modern stellar evolution calculations, we show just an example of how their synergy might explain some puzzling observations of chemical abundances in star clusters. Figure 33 displays the trend of the surface Li abundance as a function of T eff for a sample of MS stars in the approximately 600 Myr old Hyades open cluster [225,226]. It is easy to notice the so-called Li-dip around 6600 K, a sharp local minimum of the Li abundance, which cannot be explained by standard stellar evolution models including only convective mixing (the Li-dip is observed also in other open clusters of different ages, at similar temperature). Lithium is a very fragile element that is burned at temperatures of approximately 2.5 × 10 6 K. In solar metallicity MS models calculated with just convective mixing, and T eff between approximately 6400 and 7500 K, the Li burning temperature is well below the base of the convective envelope, hence the Li abundance is constant down to the radiative layers where the temperature reaches approximately 2.5 × 10 6 K. This implies that these models-which also are not expected to experience pre-MS Li depletion-predict a constant surface Li with T eff in this temperature range, contrary to observations.
To explain these observations of Li abundances, we outline below the scenario proposed by Talon & Charbonnel [227], which involves the combination of several of the element transport mechanisms discussed above. 28 To this purpose, we have divided the temperature range covered by the data of figure 33 into four contiguous ranges: (i) Region 1. These stars have shallow convective envelopes, not efficient at generating a surface magnetic field via a dynamo process, hence the surface is not slowed down by magnetic winds.  [225,226]. Inverted open triangles denote upper limits. The four different regions marked in the diagram are discussed in the text.
As a result, rotational mixing is expected to be just about sufficient to counteract atomic diffusion of Li below the convective envelope (e.g. counteracts the creation of a Li gradient right below the convective boundary where Li diffused from the convective region accumulates). (ii) Region 2. Moving towards lower temperatures, the convective envelopes deepen, and dynamogenerated weak surface magnetic fields are expected now to start slowing down the outer layers. The increased shear strenghtens the rotational mixing and Li depletion increases with decreasing T eff . This happens because the more vigorous mixing is trying now to erase the gradient between the Li-burning region (no Li) and Li-rich convective layers. (iii) Region 3. Stars on the cool side of the Li-dip have even deeper convective envelopes, that sustain a very efficient dynamo and slow down even more the external layers. At the same time, these convective layers are now very efficient at generating IGWs that redistribute momentum (driving the star towards solid body rotation) and reduce the efficiency of rotational mixing, inducing an increase in surface Li with decreasing T eff . Calculations based on the approximations described before [227] show that the expected efficiency of IGW-induced angular momentum transport has the required dependence on stellar mass (T eff ). (iv) Region 4. At these low T eff , the convective envelope is deep enough to reach Li burning temperatures, causing an increasingly larger depletion with decreasing T eff .
Measurements of projected rotation velocities in this cluster display a decreasing average velocity with decreasing T eff , consistent with an increased efficiency of magnetic braking when moving towards lower stellar masses [229].

Conclusion
It is clear that the description of complex physical processes like turbulent convection, semiconvection, thermohaline mixing and rotation, implemented in stellar evolution calculations is necessarily simplified, with the predictive power in some cases hampered by the use of several free parameters of uncertain calibration. Also, the effect of interactions among the various instabilities in rotating stars-which usually are considered as independent-has to be fully explored [230].
It is also fair to say that, in many cases, the approaches used in stellar evolution models are probably reaching their limits, and further developments of multidimensional hydrodynamical simulations are crucial to progress in the field. As we have discussed, there are already some constraints that current hydrodynamical simulations pose to the element transport processes efficient in stars. Even though there is clearly still a long way to go to complete stellar interior modelling with numerical hydrodynamics, physical insights provided by computer simulations are invaluable in improving our understanding of element transport processes [3,231]. Just to give two examples, recent results from hydrodynamical simulations provide new indications about the way to go to replace the MLT in stellar modelling [232], as well as how to implement a more consistent description of overshooting in different regimes [40].
At the same time, the booming field of asteroseismology is starting to provide new information on the efficiency of element transport processes, that can at the very least be used to add further constraints to our current recipes. Two perfect examples are recent works on mixing in low-mass core He-burning stars [233,234], and angular momentum transport in RGB stars [235,236].
The hope is that a synergy between clues from hydrodynamical simulations and the powerful constraints coming from asteroseismic analyses will help to improve our description of element transport in stellar model computations and improve their predicting power.