{'fraction_non_alphanumeric': 0.03482452027866196, 'fraction_numerical': 0.013020126500028787, 'mean_word_length': 4.617872655022641, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Double-diffusive processes in stellar astrophysics

The past 20 years have witnessed a renewal of interest in the subject of double-diffusive processes in astrophysics, and their impact on stellar evolution. This lecture aims to summarize the state of the field as of early 2019, although the reader should bear in mind that it is rapidly evolving. An Annual Review of Fluid Mechanics article entitled"Double-diffusive convection at low Prandtl number"(Garaud, 2018) contains a reasonably comprehensive review of the topic, up to the summer of 2017. I focus here on presenting what I hope are clear derivations of some of the most important results with an astrophysical audience in mind, and discuss their implications for stellar evolution, both in an observational context, and in relation to previous work on the subject.


Introduction
Double-diffusive instabilities were first discovered in the context of physical oceanography by a group of scientists from the Woods Hole Oceanographic Institution [2,3]. Stern in particular realized that a region of the ocean with a stable density stratification can nevertheless undergo fluid instabilities because the water density depends on both temperature and salinity, which diffuse at different rates. To see why that may be the case, consider first a scenario in which temperature is stably stratified (with temperature increasing upward) and salinity is unstably stratified (with salinity increasing upward), as shown in Figure 1a. This could correspond for instance to the near-surface stratification of the tropical ocean. Overall, the density decreases upward because the stabilizing temperature stratification is stronger than the destabilizing salt stratification. As such, this fluid is stable to standard convection: a parcel of fluid forcibly moved downward is less dense than its surroundings, and in the absence of any diffusion, would experience a buoyancy force that pushes it back upwards. However, if the displaced parcel is small enough, its temperature rapidly adjusts to the surroundings, but its salinity does not (because salt diffuses 100 times slower than temperature). As a result, the small parcel becomes denser than its surroundings (because it has the same temperature but is saltier) and continues to sink. A similar argument can be made for small fluid parcels initially displaced upward, that continue to rise.
This instability, called the salt fingering instability because of its tendency to produce long thin fingers of salty or fresh water, causes a net vertical transport of salt and heat downward. It is one of a few different kinds of thermohaline (or more generally thermocompositional) double-diffusive instabilities that have been discovered since the late 1950s. These include the oscillatory double-diffusive instability (ODDC), which occurs when the signs of the gradients Figure 1. Pictorial description of the two main double-diffusive instability mechanisms. The color represents potential temperature (or entropy) and the symbol density represents the concentration of a denser species. Left: the fingering instability. A parcel of fluid displaced downward rapidly equilibrates in temperature, but retains its high concentration. As a result it is denser than the surroundings and continues to sink. Right: the ODDC instability. In the absence of any background thermal stratification, a displaced parcel would just oscillate up and down (dashed line). Thermal diffusion at extreme phases of oscillation cools down or warms up the parcel, which amplifies the oscillation (solid line).
of temperature and salinity are both reversed [3,4], and the intrusive instability, which occurs when horizontal gradients of temperature and salinity are also present [5].
The idea of double-diffusive instabilities and double-diffusive mixing soon spread to the astrophysical community, thanks in part to the role of the Geophysical Fluid Dynamics summer program at the Woods Hole Oceanographic Institution, that was attended in the 1960s by astrophysicists such as Edward Spiegel, Shoji Kato, Jean-Paul Zahn, Douglas Gough, and others. In stellar interiors, the role of salt is replaced by any chemical species that has a higher mean molecular weight than its environment, as for instance He in comparison with H, or any heavier species (especially C, N, O, Si, Ni, Fe, etc. ) within a H-He mixture. In addition, whether the fluid is thermally stably stratified or not depends on the sign of the gradient of potential temperature (or equivalently, the sign of the entropy gradient) rather than on the sign of the temperature gradient alone. This accounts for adiabatic cooling or heating as the fluid parcel expands or shrinks to adapt to the local hydrostatic pressure. Fingering instabilities can take place in regions that are thermally stably stratified (i.e. radiative zones) with unstable compositional gradients. The process is usually referred to as thermohaline convection in stars. ODDC takes place in thermally unstable regions (as established by the Schwarzschild criterion) that have a sufficiently large compositional gradient to be Ledoux-stable. These regions are commonly called semiconvective, and their existence was first discussed in the late 1950s [6]. The relationship between semiconvection and ODDC was only later clarified by Kato [7] and Spiegel [8].
The physics of the ODDC instability, which takes place in these semiconvective regions, are summarized in Figure 1b. If, as a thought-experiment, we ignore the unstable entropy stratification and only consider adiabatic perturbations, the fluid is stably stratified by the compositional gradient and supports internal gravity waves. As such, a displaced parcel would just bob up and down without change of amplitude owing to its own restoring buoyancy force. The temperature stratification can destabilize this oscillation, however. Thermal diffusion causes a sufficiently small parcel to adjust to the local background temperature, which then warms it up (lowering its density) when the parcel is low, and cools it down (increasing its density) when it is high. This effect enhances the buoyancy force, and causes a gradual amplification of the oscillation, thus driving the instability.
Both fingering convection and ODDC/semiconvection were popularized in stellar astrophysics in the 1970s and 1980s. Compositional transport models were proposed for both instabilities [9][10][11][12][13] and many of them are still widely used in stellar evolution codes today. Until this last decade, however, none of these transport models had ever been tested, because laboratory experiments and numerical experiments at relevant parameters are prohibitively difficult to perform. As such, I would argue that any result ever obtained in the field of stellar astrophysics that strongly relies on double-diffusive mixing (including in particular fingering convection and semiconvection) should be taken with a very healthy dose of skepticism, and ought to be revisited in the light of the recent numerical and theoretical developments I will now proceed to review.
In what follows, I will first lay out for completeness and pedagogical purposes the simplest possible model setup in which to study double-diffusive instabilities, and will review its stability properties (Section 2). I will then describe in more depth recent results on mixing by the fingering instability (Section 3) and the ODDC instability (Section 4), and discuss their implications. I conclude with a list of model caveats, and a summary of the more recent results obtained including the effects of shear, rotation, and magnetic fields in Section 5.

Linear instability of double-diffusive systems 2.1 Model setup
In Section 1, I briefly described how the fingering and ODDC instabilities proceed, and argued that both cases involved the need for small parcels of fluids so thermal diffusion can equalize the temperature within the parcel with that of the background. The fact that doublediffusive instabilities have to be small-scale can therefore be seen as a defining property of double-diffusive convection. The lengthscale of basic double-diffusive fluid motions in both fingering and ODDC situations can usually be estimated as being of order 10d [3,14], where where κ T and ν are the thermal diffusivity and kinematic viscosity of the fluid, and N 2 T = −δg(∇ T − ∇ ad )/H p is the square of the Brunt-Väisälä frequency associated with the entropy stratification only (i.e. ignoring compositional stratification). All the quantities used here are defined with standard notations in stellar astrophysics [15,16]. Note that N 2 T can be negative if the system is unstably stratified (∇ T − ∇ ad > 0). While κ T , ν and N 2 T can vary a lot from star to star, and between the center and surface of a given star, the 1/4 power implies that the double-diffusive scale 10d itself does not vary much, taking values around 100 meters in non-degenerate regions of stars, down to about 1 centimeter in degenerate regions (assuming the fingering region extends there). We see that in all limits, d is much smaller than the stellar radius.
The fact that these are small-scale diffusively driven instabilities is both a blessing and a curse from a numerical standpoint. It implies that all diffusive lengthscales must be fully resolved at all times, and that any attempt to resolve them in global hydrodynamical models of stars is futile. Double-diffusive instabilities should also never be studied with Large-Eddy Simulations or any other numerical technique involving subgrid scale parametrizations, because they rely on microscopic diffusivities to exist. On the plus side, this also implies that these instabilities are best modeled in small computational domains that can ignore complex stellar physics. In particular, the effects of curvature, compressibility and the nonlinearity of the equation of state can almost always be neglected for such small-scale instabilities. Instead, they can adequately be studied with the Boussinesq approximation for weakly compressible gases [17].  Figure 2. Illustration of the model setup. A local Cartesian domain is modeled, near radius r m , with the vertical defining the z-axis. In the absence of rotation the orientation of x and y are arbitrary. If the system is rotating, x is aligned in the azimuthal direction.
To do so, let's consider from now on a small region of a star, located around a radius r m , with mean temperature, density and pressure T m , ρ m and p m . This region is modeled in a local Cartesian reference frame with coordinates (x, y, z) (see Figure 2), where the local gravity defines the vertical axis: g = −ge z . As such, we can define z = r − r m , for instance. Following the Boussinesq assumptions, the computational domain size must be much smaller than any of the local pressure, density, and temperature scaleheights. This guarantees that the background temperature profile T 0 (z) can be linearized close to r m , so that T 0 (z) = T m + zT 0z where T 0z is the locally constant temperature gradient 1 . Similarly, we assume that the background compositional field profile C 0 (z) (which could for instance taken to be the mean molecular weight, or the mass fraction of a particular chemical species), can be written as C 0 (z) = C m + zC 0z where C m is the mean value in the domain, and C 0z is the local gradient. Fluid motions, written as u = (u, v, w), as well as fluctuations of temperature T , density ρ and composition C, evolve around that background state following the Spiegel-Veronis-Boussinesq equations [17], namely ρ m ∂u ∂t + u · ∇u = −∇p − ρge z + ρ m ν∇ 2 u, ∂T ∂t + u · ∇T + w(T 0z − T ad,z ) = κ T ∇ 2 T, ∂C ∂t + u · ∇C + wC 0z = κ C ∇ 2 C, (4) ρ ρ m = −αT + βC and ∇ · u = 0, where κ C is the compositional diffusivity, and the coefficients of the linearized equation of state are α = −ρ −1 m (∂ρ/∂T ) C m ,p m and β = ρ −1 m (∂ρ/∂C) T m ,p m taken at r m . We see that the temperature gradient does not appear alone, but is instead combined with the adiabatic temperature gradient T ad,z = −g/c p (where c p is the specific heat at constant pressure) to form the potential temperature gradient (T 0z − T ad,z ), which is the relevant one in stellar interiors as discussed earlier. With this model setup, we can conveniently require all perturbations u, T , p and C to be triply periodic in the computational domain, which enables us to study the dynamics of double-diffusive instabilities far from any solid boundaries (that would be unphysical in a star). The set of equations (2)-(5) is commonly used to study double-diffusive con-vection with given background temperature and composition gradients far from boundaries [4,18,19].
As an illustrative side-note, these equations can very easily be used as is to look into the stability of standard (non-diffusive) multicomponent convection, and recover the well-known Ledoux criterion for instability. This is a useful warm-up exercise for the more complicated double-diffusive case. Let's first linearize the equations, substitute the equation of state into the momentum equation, and ignore all diffusive terms. We obtain (in addition to the incompressibility condition) ∂T ∂t + w(T 0z − T ad,z ) = 0 and ∂C ∂t + wC 0z = 0.
Assuming an ansatz of the form q(x, y, z, t) =q exp(ilx + imy + ikz + λt) for each of the variables, yields Assuming without loss of generality that l 0 (otherwise switch the x and y variables), we can eliminate pressure from the first two equations to getṽ = (m/l)ũ, and so from incompressibilityũ = −lkw/(l 2 + m 2 ). We can then eliminate pressure between theũ andw equation, and substituteũ,T andC in terms ofw in the result. Once completed, the process results in an equation for λ only, because (as required from linear theory) the amplitudew cancels out of the equation. We then have where N 2 = N 2 T + N 2 C and N 2 C = −βgC 0z is the contribution to the Brunt-Väisälä frequency due to the compositional stratification. With the form of the ansatz selected above, instability occurs when solutions of this equation exist for which the real part of λ is strictly positive. This is the case when αg(T 0z − T ad,z ) − βgC 0z = N 2 T + N 2 C < 0, which is equivalent to the Ledoux criterion as required.

Linear stability properties of homogeneous double-diffusive systems
When all the diffusion terms are taken into account to model double-diffusive instabilities, it is significantly more convenient to non-dimensionalize the equations first. Following standard practices in the field (see [19] for instance), the unit lengthscale from here on is d, the unit timescale is the thermal diffusion timescale across d, namely d 2 /κ T , the unit velocity is κ T /d, the unit temperature is d|T 0z − T ad,z | and the unit composition is αd|T 0z − T ad,z |/β. Note that we use |T 0z − T ad,z | instead of T 0z − T ad,z to ensure that all units are positive. Letting, e.g.
where the + sign applies in (12) and (13) for fingering convection, and the − sign applies for ODDC. It is rather remarkable to see that the only mathematical difference between these two instabilities is the sign difference in theŵ term in these equations. Since the notation is somewhat heavy as is, in what follows the hats on the non-dimensional independent variables x, y, z, t and on the operator ∇ will be dropped, but those on the non-dimensional dependent variables will be kept for clarity. Three nondimensional parameters have appeared in these equations, namely the Prandtl number Pr = ν/κ T , the diffusivity ratio τ = κ C /κ T , and the density ratio where the last expression has been written assuming C is the mean molecular µ, and using standard definitions in stellar astrophysics [15,16] for δ, φ and ∇ T − ∇ ad and ∇ µ . The Prandtl number and diffusivity ratio depend principally on the nature of the fluid considered, and their dependence on thermodynamical quantities such as temperature, pressure, composition, etc. is neglected as part of the Boussinesq approximation. In salt water for . Note that these are just indicative order of magnitudes, and the user should compute these numbers more accurately for the application of their choice. Examples of such computations are given in [20] for WDs, MS stars and RGB stars.
The density ratio on the other hand measures the relative strengths of the potential temperature and compositional stratifications. With this definition, a Ledoux-neutral stratification has R 0 = 1. As we will now demonstrate, the density ratio is the most important parameter controlling the dynamics of both fingering convection and ODDC.
The linear stability of a doubly-stratified system satisfying equations (11)-(13) can be computed more-or-less as we did before for the non-diffusive case. We first assume a similar ansatz (bearing in mind everything is now non-dimensional),q(x, y, z, t) =q exp(ilx + imy + ikz +λt). Substituting that ansatz into the linearized equations, and successively eliminating p,ṽ,ũ,T ,C and finallyw, results in the following cubic for the growth rateλ [4]: wherek = (l,m,k) is the total wavevector,k h = l2 +m 2 is the horizontal wavenumber, and where the + sign applies for fingering convection, and the − sign applies for ODDC. As in the non-diffusive case, the existence of solutions of this cubic with positive real part demonstrates instability.
Despite the strong similarity of the cubics obtained in the fingering and ODDC cases, respectively, the solutions can be quite different. Using standard properties of cubics, it can be shown that fastest-growing modes in the fingering regime (with R 0 > 1) are real while they are complex in the ODDC regime (with R 0 < 1). For this reason, we now discuss the two cases separately.

Linear stability properties of the fingering regime
To study the fingering regime, we consider (15) with the + sign. Solutions are real so the condition for marginal stability can simply be written asλ = 0, which implies that the critical density ratio for instability for a given mode with wavenumberk is The largest possible R 0,c (k) is obtained whenk = 0, and is equal to R f c = τ −1 . We have therefore demonstrated that the fingering instability exists when 1 ≤ R 0 ≤ τ −1 . The system is unstable to standard overturning according to the Ledoux criterion when R 0 < 1, and stable when R 0 > τ −1 . These findings are not too surprising: in this regime, R 0 can be interpreted as the ratio of the stabilizing stratification to the destabilizing one, so the larger R 0 is, the more stable the system is. We also see that the range of instability depends only on the value of the diffusivity ratio τ = κ C /κ T . If temperature and composition diffuse at the same rate, then fingering is not possible. On the other hand, the smaller τ is, the larger the range of density ratios for which instability can exist. Because τ is usually asymptotically small in stars, this implies that even a very small inverse µ gradient can destabilize a radiative zone. With a little work, it can be shown [19] that the modes withk = 0 are always the most rapidly-growing ones, and are often called elevator modes for obvious reasons: the flow within these elevator modes is strictly vertical, and all the quantities of interest are invariant with z. In the limit of low Prandtl number and low diffusivity ratio appropriate for stellar interiors, and close to marginal stability (so R 0 → τ −1 and bothλ and |k| are small), we can roughly estimateλ for these elevator modes by neglecting the cubic and quadratic terms in (15) (see [9]):λ In reality, however, the horizontal wavenumberk h of the fastest growing elevator mode also varies with R 0 close to marginal stability, so this expression does not directly tell us what λ is for that mode. To do so, one needs to perform a formal asymptotic expansion of (15) in the limit of low Pr and τ. This was done by Brown et al. [21], who derived a number of approximate analytical solutions for the growth rateλ of the fastest-growing modes as a function of R 0 , Pr and τ, in the stellar limit (Pr, τ 1). The reader is referred to their Appendix B for details. Among other results, they find that for R 0 not too close to 1 nor too close to marginal stability (i.e. 1 R 0 τ −1 , which seems to be a reasonable limit for stars),λ while the horizontal wavenumber of the fastest-growing mode is always of order unity, witĥ Dimensionally, this implies that the wavelength of these fingers is of order 2πd (corresponding to one up-flowing and one down-flowing finger), where d was estimated in (1), and that their growth rate is which is usually substantially smaller than the buoyancy frequency unless R 0 is close to 1.

Linear stability properties of the ODDC regime
To determine the stability properties of ODDC, we study the solutions of (15) with the − sign.
In that case, it is common to use the inverse density ratio R −1 0 as the relevant parameter; R −1 0 is now the ratio of the stabilizing compositional stratification to the destabilizing potential temperature stratification, with R −1 0 < 1 being Ledoux-unstable (hence convective in the standard sense), while R −1 0 > 1 is Ledoux-stable. Increasing R −1 0 correspond to increasingly stable systems. The criterion for marginal stability to ODDC is a little more difficult to derive than in the fingering case, becauseλ is not real. So we letλ =λ R + iλ I and first rewrite (15) in terms ofλ R only. The same argument as in the fingering case can be used to show that the elevator modes are still the fastest-growing modes, so focussing on these only, we letk = 0. With a little work (see the Appendix A.1 in [22]), we obtain Settingλ R = 0, we then find that the critical inverse density ratio for instability is achieved fork h → 0 (as in the fingering case). This shows that ODDC can only occur for Since Pr τ in geophysical environments, (Pr + 1)/(Pr + τ) 1. As a result, the range of R −1 0 for which a doubly stratified system on Earth may be linearly unstable to ODDC is very small and almost never naturally realized. By contrast, both Pr and τ are very small in stars, so the range of instability to ODDC can be very substantial.
Applying the same arguments as in the fingering case, we can roughly estimateλ R close to marginal stability (i.e. when R −1 0 → (Pr + 1)/(Pr + τ) where bothλ R andk h are small) by neglecting the quadratic and cubic terms in (21). In the stellar limit where Pr, τ 1, we then get [12],λ for sufficiently large R −1 0 (but not too close to marginal stability). Again, this should be viewed as a rough estimate, and does not provide information on the wavenumberk h for the fastest-growing modes. A more formal asymptotic analysis needs to be carried out to get this information. Details can be found in Appendix A.2 of Mirouh et al. [22]. The solutions have an analytical but non-trivial dependence on R −1 0 , whose expression is not particularly illuminating. At best, one can state thatλ R is proportional 2 to Pr, so from a dimensional point of view, we find that λ R ∝ √ PrN T . Finally, it is also relatively easy to show that for Pr, τ 1, the imaginary part ofλ is well-approximated by the local buoyancy frequency including the compositional stratification, i.e.

Summary for linear stability of double-diffusive systems
An illustrative summary of the findings from linear stability analysis is presented in Figure 3.
In the absence of compositional stratification, the threshold between stability and instability to overturning convection is simply given by the Schwarzschild criterion, namely R 0 = 0 or equivalently ∇ T = ∇ ad . With compositional stratification, instability to overturning convection is set by the Ledoux criterion, R 0 = 1 or equivalently δ(∇ T − ∇ ad ) = φ∇ µ . We see that double-diffusive effects destabilize regions that are Ledoux-stable. In the presence of a destabilizing compositional gradient, fingering instabilities are excited in a wide region of parameter space that would normally be purely radiative. In the presence of a stabilizing compositional gradient, ODDC is excited almost everywhere in the region of parameter space bound by the Schwarzschild criterion on one side, and the Ledoux criterion on the other. In all cases, the fastest-growing modes of instability are elevator modes, i.e. vertically-invariant structures, whose horizontal size is O(10d).
Overturning convection Stable system

Overturning convection
Stable system While this analysis has been performed for an idealized system without any added physics, we are of course also interested in knowing whether rotation, shear, magnetic fields, and other dynamics, could affect these results. Interestingly, neither rotation, magnetic fields, nor shear seem to affect the range of instability of fingering convection [23][24][25], because there always appears to be a way for the elevator modes to develop. Horizontal gradients, however, can extend the unstable fingering range significantly [5,26]. Gravitational settling may induce as similar effect if the settling velocity of individual species approaches the value κ T /d [27,28], though this has not yet been studied in the astrophysical context. To my knowledge, there is no evidence for subcritical instabilities (i.e. instabilities that only develop given specific finite-amplitude initial conditions instead of from infinitesimal perturbations) in fingering-prone stratified fluids.
The impact of added dynamics on ODDC has not yet been studied as exhaustively as in the fingering case. Rotation does not affect the unstable range [29], but shear does (at least in geophysics), through the newly discovered thermo-shear instability [30]. Whether that instability is relevant at astrophysical parameters has not yet been established. The effect of magnetic fields on the linear stability of ODDC was investigated by Stevenson [31], who found that it does not affect the overall unstable range, but can change the nature of the unstable modes that are present. Finally, by contrast with the fingering case, ODDC is known to have subcritical branches of instability, at least in the geophysically-relevant region of parameter space [32]. This remains to be studied more extensively at low Prandtl number [33].
Finally, note that all of these results have been stated in the context of the model setup described in Section 2.1. The presence of physical boundaries can strongly affect the instability range, especially when the size of the domain is not very large compared with the natural double-diffusive scale 10d. For bounded domains the linear stability problem needs to be investigated on a case-by-case basis, and depends on the domain size, shape, applied boundary conditions, and any added physics.

Where is fingering convection likely to occur in stars?
The linear stability analysis presented above clearly shows that fingering instabilities can take place in radiative stellar regions which are Ledoux-stable, but subject to an unstable composition gradient, even if the latter is very weak. Such gradients can arise in stellar interiors in a number of scenarios, all of which have important observational implications.
A common mechanism to drive fingering convection involves the accretion of high µ material at the surface of a star, either from infalling planets or planetary debris, or from material transferred from a more evolved companion star. Planetary infall for instance can increase the apparent metallicity of the host star; this has been invoked in turn as a possible explanation to the planet-metallicity correlation [34], and to the existence of metal-rich WDs [35,36]. It was rapidly realized however that this inverse µ-gradient would create a region below the surface of the star that is unstable to fingering convection [37], which would then significantly shorten the residence time of high-µ chemical species near the surface, with implications for derived accretion rates and/or surface abundances of light elements in both MS stars [38,39] and WDs [40][41][42]. Similar arguments have been put forward in the context of binary mass transfer, with possible implications for carbon-enhanced metal poor stars, for instance [43][44][45][46] and nova outbursts [47].
Fingering regions can also appear deep within a star from off-center nuclear burning in RGB and AGB stars. In RGB stars, 3 He burning at the outskirts of the H-burning layer can produce lower-µ material, causing that region and the regions above it to become fingeringunstable [9,48]. Instabilities associated with this inverse µ gradient have been proposed as a possible solution to the peculiar surface abundances of RGB stars beyond the luminosity bump [49,50]. Fingering is also thought to occur below the carbon-burning shell of super-AGB stars [51], and to affect the carbon flame propagation by mixing the available fuel. It may also impact the properties of hybrid C/O/Ne WDs, whose electron-to-baryon fraction becomes unstable to convection and fingering convection as the star cools [52,53].
Finally, the microscopic segregation of various chemical species by gravitational settling and radiative levitation can produce chemically peculiar layers below the surface of intermediate-mass stars, that are prone to fingering instabilities [54][55][56]. Not accounting for mixing by fingering, the accumulation of such elements can be so strong that it leads to the formation of intermediate convection zones. Fingering instabilities, however, may develop much earlier and prevent the accumulation from reaching such extreme levels [55].

Where is ODDC/semiconvection likely to occur is stars?
We also learned using linear theory that ODDC develops in regions that are traditionally called semiconvective, i.e. regions which are unstable according to the Schwarzschild criterion (∇ T − ∇ ad > 0), but stable according to the Ledoux criterion (φ∇ µ > δ(∇ T − ∇ ad )).
Semiconvective regions are often found adjacent to convective cores and are caused by the development of stabilizing µ gradients. There are two stellar mass ranges in which this occurs, and the reason for the presence of the µ gradient differs in the two cases. For reviews on the topic, see for instance [8,12,13,57,58].
In intermediate mass MS stars, with masses above 1M and below 3M , nuclear reactions can extend far beyond the radius of the convective core, owing to the relatively weak dependence of the pp-chain reaction rates on temperature. This implies that He is slowly generated outside the core, at a rate that decreases with radius. This creates a gentle µ gradient, that partially inhibits convection and causes the convective core to be smaller than what it would be according to the Schwarzschild criterion. A semiconvective region develops between this so-called Schwarzschild core radius and the Ledoux core radius, see Figure 4a.
In the core of higher mass stars (with masses greater than about 5M ) where the CNO cycle dominates, the temperature dependence of the nuclear reaction rates is so strong that the latter only takes place deep within the convective core. As such, the mechanism described above for creating µ-gradient near the core of intermediate mass stars does not work. However, the convective core shrinks with time as H is converted into He, because the opacity of the material in that region is dominated by electron scattering and is simply proportional to 1 + X, where X is the hydrogen mass fraction. With a reduced opacity, more of the energy can be transported radiatively, and the convection zone retreats. This leaves behind concentric shells of increasingly high µ material, and therefore gradually builds up a substantial µ gradient outside this core (see Figure 4b). That region can be unstable to ODDC 3 .
Finally, semiconvection zones detached from nearby convective zones are sometimes found in models of high mass stars (15M or higher), see for instance [59]. However, whether they exist in stars or not remains to be established, because their presence or absence in models is quite sensitive to the semiconvective mixing prescription used.  Figure 4. Illustration of the two mechanisms for forming semiconvective regions at the edge of a convective core. The red line represents the nuclear energy generation rate, and the other colored lines represent the µ(r) profile, with blue at early times, then green, then orange at latest time. The vertical dashed line marks the edge of the convective core according to the Ledoux criterion. Left: In intermediate mass stars, nuclear reactions extend outside the core, and their rate depends on temperature. A µ gradient thus forms in situ. Right: In high-mass stars, nuclear generation only takes place well within the core. However the core shrinks with time (see text for detail) and leaves behind material of increasing µ. A µ-gradient gradually forms as a result. In both figures r L marks the edge of the core according to the Ledoux criterion. The ODDC can exist just outside of it.

Beyond linear theory
Having established the conditions under which fingering convection and ODDC may develop in stars, we can now ask the more important quantitative question of how much mixing they cause, and how this affects stellar evolution. To do so requires connecting fluid dynamics, which governs the double-diffusive instabilities on the small scales / short timescales, and stellar modeling, whose equations govern the structure and evolution of stars on the largescales / long timescales.
In stellar evolution codes, turbulent transport by double-diffusive convection is traditionally modeled as a turbulent diffusion process. Assuming the star is spherically-symmetric, the conservation law for the evolution of the concentration C of a particular chemical species in a given mass shell is where the time-derivative is a Lagrangian derivative following the shell, ρ is its density, ρC is the total mass of the particular chemical element considered, ρF C is the mass flux, and (D(ρC)/Dt) nuc is the rate of change due to nuclear (or chemical, if relevant) reactions. The assumption that turbulent transport takes a diffusive 4 form implies that the compositional mass flux F C should be downgradient, with where κ C is the microscopic diffusivity and D C is the turbulent diffusivity, both of which have units of cm 2 /s in cgs. The assumption of spherical symmetry then implies that which can be expressed in mass coordinates m as This is the formula typically implemented in stellar evolution codes for compositional transport, with a turbulent mixing coefficient D C that depends on the instability driving the turbulence. It is worth remembering, however, that turbulent transport does not always take a diffusive form, so one should ideally first question whether that assumption is correct, before attempting to model D C . The case of turbulent heat transport by double-diffusive instabilities could in principle be treated in a similar way, but in practice this is rarely ever done, for two reasons. The first is that very few stellar evolution codes actually evolve the temperature profile in the star, preferring instead to solve for it at each timestep knowing the stellar luminosity. The second is that heat transport in fingering convection and semiconvection is usually negligible (with some notable exceptions discussed later), because temperature has to diffuse for the instability to exist in the first place.
The rest of this lecture will be therefore be dedicated to presenting models for D C for double-diffusive instabilities, with the fingering case discussed in Section 3, and the case of ODDC in Section 4. The following Table summarizes the important non-dimensional parameters that characterize the properties of double-diffusive convection. All will be extensively used in this lecture. Input parameters such as Pr, τ and R 0 have already been discussed in this Section. Output non-dimensional parameters based on D C (and other turbulent mixing processes) will be introduced and used in Sections 3 and 4, to help quantify turbulent transport of heat and composition in the system, both in comparison with the respective diffusive transport rates, or in comparison with one another. Table 1. Non-dimensional parameters governing double-diffusive mixing, written as the ratio of dimensional quantities. The dimensional turbulent fluxes of temperature and composition are F T = wT and F C = wC , where the notation · denotes a volume average. Assuming turbulent transport is diffusive implies that F C = −D C C 0z in this notation.

Name
Symbol and Formula Interpretation Input params.
Ratio of turbulent fluxes Total flux ratio

Ratio of total fluxes
3 Mixing by fingering convection

Traditional models of "thermohaline" mixing
The first turbulent diffusion model proposed for mixing by fingering convection in astrophysics was put forward by Ulrich [9], and is based on a very simple dimensional analysis. Noting, as it is common to do so in astrophysics, that the turbulent diffusion coefficient D C has the units of a velocity times a length, or equivalently, of a length squared divided by time, it is reasonable to assume that the one appropriate for fingering convection should be expressed as where λ is the dimensional growth rate of the fastest-growing fingers, which is related to the non-dimensional growth rate discussed in Section 2 via λ = κ Tλ /d 2 . Using the estimate from (17) for the growth rate of fingers close to marginal stability, we get Assuming thatk 2 h remains close to unity, then we can write In the limit where R 0 is not too close to marginal stability (which is somewhat inconsistent with the assumption made above, but let's ignore this for now) then 1 − R 0 τ 1 and we recover the formula proposed by Ulrich [9], with a proportionality constant he argues should be of order 700: A very similar expression was later obtained from rather different arguments by Kippenhahn et al. [11], who argued that We see by comparison with (31) that this expression should only be valid in the limit where 1 R 0 τ −1 . The coefficient C KRT in this model is argued to be much smaller than C U , taking a proposed value of 12.
We can already see, however, that both models as they are presented fail to account for the stabilization of the system to fingering instabilities beyond the threshold R 0 = τ −1 . As such, we expect that they should largely overestimate the true mixing coefficient as R 0 approaches and/or exceeds that threshold. A better model, that would at least be consistent with linear stability theory, is given in (31) and was (to my knowledge) first derived by Denissenkov [60] (see his equation 15). Written in terms of the more standard astrophysical notations, the model becomes assuming the compositional field C is simply the mean molecular weight. That model, however, is ill-posed (with D fing → ∞) as R 0 tends to unity, or in other words, as we approach the Ledoux criterion for convective instability. This is obviously not physically plausible, so the model should not be used when R 0 → 1. This ill-posedness is not surprising since the model was derived in the first place assuming that R 0 is close to marginal stability, i.e. very large; better models (see below) have since been proposed to address the problem.

Numerical simulations of small-scale fingering convection
Until recently, no robust evidence for or against the adequacy of the various expressions for D fing listed above existed. As mentioned in the introduction, laboratory experiments at low Prandtl numbers are almost impossible, because so few low Prandtl number fluids exist on Earth, and those that do are very expensive and notoriously difficult to manipulate (e.g. liquid mercury, liquid lithium, liquid potassium). Numerical experiments are also challenging. Indeed, the computational domain needs to contain at least 5-10 fingers in each direction for good statistics. Furthermore, we saw that the scale of the fingers is related to the thermal diffusion scale, while the scale of boundary layers between the fingers is dictated by the viscous and compositional diffusion scales, which are asymptotically small under stellar conditions (since Pr = ν/κ T 1 and τ = κ C /κ T 1). All of these scales need to be resolved to adequately model fingering convection, which poses a hard computational challenge. Finally, it has recently been demonstrated that fingering convection at low Prandtl number cannot be modeled in 2D, as it develops unphysical pathological behavior [61]. In 2D simulations, artificial horizontal shear layers appear spontaneously from the fingering instability, and in turn affect the fingering structures. These shear layers are not present in 3D domains 5 , as long as the third dimension is wide enough (i.e at least 2 finger widths), see Figure 5. 5 at least not as ubiquitously and with such strong amplitude.  . Comparison between 2D simulations and thin-domain 3D simulations at the same parameter values (Pr = τ = 0.01, and R 0 = 5). In the 2D case, horizontal shear layers spontaneously develop and tilt the fingers, which in turn reduces the turbulent fluxes. This is artificial, however, and these shear layers are not present in 3D as long as the domain is thick enough. Figure adapted from [61].
Within the last decade, however, numerical simulations with a physical grid resolution of ∼ 300 3 have become routine, and state-of-the-art ones are now exceeding ∼ 3000 3 or more. With the smaller-sized runs, it is possible to systematically explore parameter space (i.e. cover the entire range of density ratios) for fingering convection at diffusivity ratios Pr and τ down to ∼ 0.01. The higher resolution simulations can be used to test models down to, e.g. Pr, τ ∼ 0.001, but only for a few values of the density ratio.
A first series of 3D DNSs of fingering convection at low Prandtl number was presented by Traxler et al. [62], using the PADDI code developed by S. Stellmach. This code is a pseudo-spectral code that solves equations (11)-(13) in a triply periodic domain, ensuring incompressibility is maintained using a standard projection method [63]. Traxler et al. limited their study to moderate Pr, τ, from 1/3 down to 1/30. Figure 6 shows representative snapshots of the compositional field in their simulations, for two density ratios (one low, close to the Ledoux threshold, and one large, close to the marginal stability threshold). Simulations with lower Pr and τ down to 1/300 made with the PADDI code are now also available [1,21] and show qualitatively similar features. Each simulation can easily be used to measure a turbulent compositional diffusivity. Indeed, in a triply-periodic domain, the horizontal average (marked with an overbar) of the dimensional composition evolution equation (4) is which reveals wC to be the vertical turbulent compositional flux (i.e. the transport of C due to advection by vertical fluid motions). While this quantity would generally depend on both z and t, in a statistically steady state and assuming thatC remains small (which we can verify is true in these simulations), wC is approximately constant. We can then obtain a good estimate of the vertical turbulent flux by taking both a time-average and a volume average of wC in the simulation, and define this as F C = wC . The same can be done for the temperature. Incidentally, it is easy to show [64] that F C and F T must be negative for statistically stationary fingering convection. Indeed, multiplying equation (4) by C, integrating over the volume, and using incompressibility and the periodicity of the boundary conditions to eliminate the nonlinear advection terms and the boundary terms, we get using one of Green's identities. Since C 0z > 0 in fingering convection, and since |∇C| 2 > 0, then wC must be negative if the turbulence is statistically stationary (i.e. if we can neglect the time-derivative). The same is true for wT , showing that both T and C are transported downward (which is not surprising, given the physical mechanism described in Section 1). The volume-averaged compositional flux can finally be used to derive D fing from each simulation as where the hats denote non-dimensional quantities (see Section 2.2). The numerical results can be used to test the traditional models of Ulrich and Kippenhahn et al. [9,11] described in Section 3.1. Writing both models as D fing = C U,KRT κ T R 0 −1 , with only the constant differing between them, the theory would predict that D fing (R 0 − 1)/κ T should be constant. The results are shown in Figure 7, for a wide range of simulations [1,21,62]. For ease of visualization the results are shown, not as a function of R 0 , but as a function of the reduced density ratio [62] which maps the entire fingering range into the interval [0, 1] regardless of τ, with r = 0 corresponding to the Ledoux criterion, and r = 1 corresponding to marginal stability. We can immediately make several important conclusions: that D fing (R 0 − 1)/κ T is not constant across the entire range, that it seems to depend on the ratio ν/κ C = Pr/τ (called the Schmidt number) and finally that the numerical results do not support the large value of C U 700 [9] but are more consistent (within the list of caveats listed) with the smaller value C KRT 12 [11]. As discussed earlier, the fact that D fing (R 0 −1)/κ T is not constant across the fingering range is not surprising, since D fing ∝ κ T R 0 −1 is both unphysically singular as R 0 → 1, and does not account for the stabilization of the system as R 0 → τ −1 (which is clearly visible in the data). The fact that the results depend principally on the Schmidt number is also expected, both on physical and mathematical grounds. Physically speaking [62], when κ T (κ C , ν), the role of the temperature fluctuations essentially becomes negligible, and the instability is driven by the compositional field. It is therefore not surprising to find that the dynamics solely depend on ν/κ C , rather than on any parameter that depends on κ T . From a mathematical point of view, it can be shown [65,66] that in the limit of κ T → 0 (and for sufficiently weak fingering) the non-dimensional governing equations (11)-(13) reduce to a new set of asymptotic equations that only depend on two parameters instead of the usual three, namely the Schmidt number, and the so-called Rayleigh ratio R = (R 0 τ) −1 . As such, the data dependence on the Schmidt number, rather than Pr and τ, is expected.
Based on their first set of numerical experiments, Traxler et al. [62] proposed a simple empirical formula for mixing by fingering convection, namely where the numbers 101, 3.6 and 1.1 were fitted to the data. This model is easy to use, and fits the results adequately except for very low r [21], where the system is only weakly stratified and close to the Ledoux threshold for overturning convection. In that limit (39) dramatically underestimates the mixing coefficient [21] .  Figure 7. Non dimensional ratio D fing (R 0 − 1)/κ T , against the reduced density ratio r, for a wide range of simulations. This quantity should be constant according to the traditional models of fingering convection [9,11] but is not in reality.

The Brown et al. 2013 model for small-scale fingering convection
In order to address the shortfalls of the Traxler et al. model for fingering convection, Brown et al. [21] revisited the problem and proposed a new theory for the mixing coefficient. Starting with the definition of D fing given in (37), we see that the key to creating a model for this coefficient is to estimate both the typical non-dimensional vertical velocity of the fingers,ŵ f , and their typical compositional perturbation,Ĉ f , so where K B is a constant that depends on the geometry of the fingers, and the typical correlation betweenŵ andĈ. Assuming that transport is controlled by the most rapidly growing fingers (which are elevator modes) and that these fingers are controlled by the linearized equations (at least until they nonlinearly saturate) we can relateŵ f andĈ f using the linearized version of (13), in which ∂/∂t is replaced by the finger growth rateλ, and horizontal derivatives are replaced by the horizontal wavenumberk h , sô Hence and the only problem remaining is to estimateŵ f . To do so requires specifying the mechanism by which the linear fingering instability saturates. Following Radko & Smith [67], Brown et al. [21] assumed that saturation occurs by parasitic shear instabilities that develop between up-flowing and down-flowing fingers, and that the fingers stop growing when the parasitic instability growth rateσ approachesλ. By dimensional analysis (or by solving the problem exactly, see Appendix A of Brown et al.), it can be shown thatσ = K σŵ fkh where K σ is a universal constant of order unity (whose value is irrelevant for reasons explained below) so the saturation condition reads where K w is another universal constant of order unity. We see that these two constants end up folded into a single one, which is the reason why their individual values are irrelevant.
The relationship betweenŵ f at saturation of the fingering instability andλ/k h was recently verified [23] against a re-analysis of the existing numerical data, and revealed that K w /K σ 2π. With this, we finally have the prediction that This can easily be tested against available data, and the combination of constants K B (K w /K σ ) 2 is found to be around 49 (implying that K B 1.24). The comparison between the model and the data is very good (though not perfect) for all cases where ν/κ C ≤ 1, which is indeed the case in stellar interiors, for all r, Pr and τ tested. For decreasing Pr and τ the good fit deteriorates somewhat, though remains reasonable (within a factor of order unity), see Figure  8.
In summary, to compute D fing one should first find the growth rate and wavenumber of the fastest-growing modes of the fingering instability at the desired parameter values (R 0 , Pr, τ), which can be done numerically using a cubic root-finding algorithm. Onceλ andk h are known, they are used in expression (44) to predict D fing . This method has already been implemented for instance in MESA. Alternatively, analytical approximations forλ andk h (found in Appendix B of Brown et al. [21]) can be used instead to speed up the process. Routines providing such estimates are also available in MESA.

Large-scale instabilities?
While the process of small-scale fingering convection in stellar environment is now arguably well-understood, this may not be the full story. Indeed, fingering-unstable regions of the ocean on Earth are often associated with thermohaline staircases, which are horizontallyinvariant stepped structures in the vertical profiles of temperature and salinity (see [19] for a review). These staircases are formed of stacked layers and interfaces. Within a layer, the density is very slightly unstably stratified and subject to large-scale convective overturning. Both temperature and salinity are almost constant within each layer as a result of the strong mixing.
In between the layers are stably stratified interfaces undergoing fingering convection. While the layers can be tens or even hundreds of meters deep, the interfaces are much shallower (tens of centimeters), so the temperature and salinity gradients across the interface are very large (hence the characteristic appearance of the staircase). When thermohaline staircases are present, vertical mixing can be enhanced by two orders of magnitude [68] compared to a similar overall stratification without staircases. It is therefore crucial to understand why and how such staircases form in the ocean, and whether similar processes may be taking place in stars.
Thankfully, significant progress in modeling the formation of oceanic thermohaline staircases (and the emergence of other large-scale dynamics) from fingering convection has been made in the past 20 years. Given the wide separation between the finger scale and the scale of the staircases, mean field hydrodynamics turns out to be a fruitful approach to the problem. In this type of theory, the large scales (such as the layer scale) are modeled exactly, while the effect of small scales (i.e. the basic fingering convection) is parameterized using some form of turbulence closure. As we saw earlier, for a given fluid (i.e. given Pr and τ), the turbulent fluxes in fingering convection only appear to depend on the local density ratio. With this in mind, it is easy to see how large-scale instabilities might develop (see Figure 9). Indeed, any large-scale perturbation in the temperature and compositional fields causes largescale modulations in the local density ratio (which depends on their gradients). This in turn modulates the temperature and compositional fluxes due to fingering, and the convergence or divergence of these fluxes can in some cases enhance the original perturbation, and in some cases suppress it. Positive feedback loops, when they exist, thus drive the growth of largescale instabilities. As we have discovered [20,63], several distinct types of positive feedback loops can in fact exist, each leading to the amplification of different kinds of perturbations.

Perturbations
in local density ratio   Figure adapted from [1]. The background image is a snapshot from a simulation exhibiting the spontaneous emergence of large-scale gravity waves from [61] (see also [20]).
To model this idea mathematically and determine whether thermocompositional layers may form in stars, we essentially follow the theory of Radko [69], and expand it to account for the diffusive contribution to the fluxes [62] (which, as we shall demonstrate, is a crucial element of the stellar problem). Since layers are horizontally invariant, we begin by taking the horizontal average of the non-dimensional equations for temperature and composition, namely ∂T ∂t where we expectw = 0 otherwise mass is not conserved. Note that the hats have been dropped from the horizontally averaged quantities to avoid crowding the notations, but everything in this section is implicitly non-dimensional. If we then define the total fluxeŝ as the sum of the turbulent flux plus the diffusive flux of each quantity, then these equations simply become the conservation laws We also define two non-dimensional quantities: the Nusselt number Nu T , and the flux ratio γ tot , as which are similar to those defined in Table 1, but expressed as the ratio of non-dimensional quantities, and including large-scale perturbationsT andC. The Nusselt number is the ratio of the total temperature flux to the diffusive (potential) temperature flux, the latter being simply equal to −(1 + ∂T /∂z) in these units. The key assumption made by the Radko [69] is that Nu T and γ tot can only depend on other local non-dimensional properties of the fluid. Aside from Pr and τ, which are fixed once the fluid is specified, the only other relevant non-dimensional quantity is the local density ratio. The latter is given by since 1 is the non-dimensional background potential temperature gradient, and R −1 0 is the nondimensional background compositional gradient. To determine the evolution of large-scale horizontally-invariant perturbations, we therefore simply evolve the equations in (47) together with equation (48), where the functions Nu T (R; Pr, τ) and γ tot (R; Pr, τ) are assumed to be known (they can be derived for instance from the Brown et al. model or from experimental data), and R given in (50).
A trivial solution of these equations exists: whenT ≡C ≡ 0, then R = R 0 , Nu T = Nu T (R 0 ; Pr, τ) and γ tot = γ tot (R 0 ; Pr, τ) are all constant, so the fluxesF T,tot andF C,tot are also constant. This defines a turbulent state that is spatially homogeneous and statistically steady -this is the basic state of small-scale fingering convection.
WhenT andC vary with z, solutions cannot be found analytically in general because the functional form of Nu T and γ tot can be quite complex, and (50) is nonlinear. Instead, we proceed by linearizing the large-scale equations around the homogeneous fingering solution described above to study its stability to the development of large-scale perturbations. Let us therefore assume thatT andC are small, so that ∂T /∂z 1 and ∂C/∂z R −1 0 . Then We continue by linearizing Nu T (R) (at fixed Pr and τ) in the vicinity of R 0 , as Using this with (47) and (48), we get after successive simplifications (from linearization) where where A γ ≡ R 0 , and where we used the fact that in the absence of any large-scale perturbations,F T,tot (R 0 ) = −Nu T (R 0 ). Finally, assuming thatT ∝ exp(iKz +Λt) and similarly forC, we can substitute these solutions into (52) and (53), and with a little algebra, obtain a quadratic equation for the non-dimensional growth rateΛ of these horizontally-invariant perturbations: where for simplicity of notation I have used Nu 0 = Nu T (R 0 ) and γ 0 = γ tot (R 0 ). Solutions for Λ that have a positive real part denote instability. When unstable, perturbations of the form T ∝ cos(Kz) or sin(Kz) (and similarly forC) grow exponentially with time, until the density profile itself has regions that are dynamically unstable (i.e. where density increases upwards). As soon as this happens, overturning convection sets in, and stacked convective layers appear separated by sharp interfaces (see Figure 10). From the form of (54), we can immediately deduce a few important properties of the solutions. First, for unstable modes to exist (i.e. solutions with positive realΛ), it is sufficient to require that A γ > 0, or equivalently, that γ tot be a decreasing function of R. This is called the γ-instability criterion, first derived by Radko [69].
Second, we see that there is no stabilization of the system for largeK. Instead, it is easy to show thatΛ is always proportional toK 2 , so the smallest scale modes grow the most rapidly, which is not particularly physical. At the root of this so-called ultraviolet catastrophe, with ρ 0 (z) + ρ(z) z Figure 10. Eigenmodes of the layering instability are sinusoidal perturbations in the density profilē ρ(z). They grow exponentially, and eventually cause the total density profile ρ 0 (z) +ρ(z) to increase with height at specific positions in the domain. Being unstable to overturning convection these regions are then fully mixed and become the layers of the staircase. Credit: S. Stellmach.
Λ ∝K 2 , is the anti-diffusive nature of the fingering fluxes when γ tot decreases with R. To see this, let us construct the evolution equation for the horizontally-averaged density perturbation (using as before a linearization of the equations around R = R 0 ) where I have linearized R to go from the second to the third line, and again usedF T,tot (R 0 ) = −Nu 0 in the last line. Since Nu 0 is always positive in fingering systems, we see that the first term behaves anti-diffusively if A γ > 0 and diffusively if A γ < 0, which confirms the statement made above. Note that in reality the mean-field equations stop being valid on scales approaching the finger scale, and should not be used in that limit. As such, this ultraviolet catastrophe is only a feature of the mean field equations but would not occur in a real system. The γ-instability model has been quantitatively validated in the geophysical context for fingering convection at high Prandtl number in both 2D [69] and 3D [70], confirming its ability to predict not only when layers should form, but also at what rate (as long as the wavenumber of the layering mode is not too large). As a result, we can use it to determine whether layers are expected to form in stars or not. To do so, we simply need to compute the function γ tot (R) and apply the γ-instability criterion.
Nondimensionally, in a homogeneous fingering simulation at density ratio R 0 , so we can use the numerically-determined fluxes ŵT and ŵĈ from e.g. Traxler et al. and Brown et al. [21,62] to compute γ tot (R 0 ) at moderate vales of Pr and τ. The results are shown in Figure 11. In all cases, we see that γ tot is a strictly increasing function of the density ratio. The same can be shown to be true at stellar values of Pr and τ, with estimates for ŵĈ and ŵT obtained using the Brown et al. model. To understand why this is the case, note that the turbulent fluxes decrease significantly when Pr and τ do, and become small compared with the diffusive fluxes [21]. As a result, for sufficiently low Pr and τ, we have γ tot R 0 τ −1 which clearly increases with density ratio. This effectively demonstrates that layering cannot be produced from the γ−instability in the stellar regime. Curiously, Brown et al. [21] did report that one of their simulations at very low density ratio spontaneously developed layers. However, these layers were not formed by the γ−instability, but instead, by the nonlinear development of large-scale internal gravity waves, that were themselves excited by the small-scale fingering. An example of such gravity waves can be seen for instance in the background snapshot of Figure 9. The spontaneous emergence of internal waves from fingering convection is another well-known large-scale instability, sometimes called the collective instability, that can also be modeled using mean field hydrodynamics [20,63,71]. The algebra associated with that calculation is however much more involved than the one outlined above to model the layering instability, so I will not derive it here. The reader is referred to the work of Traxler et al. [63] for details of the calculation (see also [20]).
Using the turbulent flux laws from Brown et al. [21] to close the mean field equations, we were able to determine that fingering convection is expected to excite internal gravity waves down to Pr, τ ∼ 10 −3 but not lower [1,20]. Furthermore, they can only develop when the density ratio is close to one (i.e. close to the Ledoux limit). This finding rules out the excitation of gravity waves from basic fingering convection in non-degenerate regions of stellar interiors (where Pr, τ ∼ 10 −6 ), but not in degenerate regions of e.g. WDs and evolved stars. Whether a fingering region would extend all the way into the degenerate region remains to be determined, however. Furthermore, even if gravity waves are indeed excited by this mechanism, they would remain relatively small scale, so would probably not be observable anyway [20]. As such, this effect seems to be more of an interesting curiosity rather than something that is likely to impact stellar evolution and observations 6 .

Conclusions for now
Having established that fingering convection is not likely to spontaneously excite larger-scale dynamics (layers or internal gravity waves) at low Prandtl number and diffusivity ratio, we can conclude that in the absence of any other dynamical process (rotation, magnetic fields, shear, etc.), fingering-induced mixing in stars remains small scale, and is well-described by the model of Brown et al. [21] (see Section 3.3). The role of rotation, magnetic fields, and shear on fingering-induced mixing is the subject of ongoing work and will be briefly discussed in Section 5.

Applications to stellar astrophysics
In Section 2.3, we saw a number of examples where fingering convection may occur in stars. Armed with a better quantitative understanding of the process, we are now equipped to answer the question of how much mixing it causes and what its observable effect on stellar evolution may be. Here, it is almost impossible to give a comprehensive review of the topic. Instead, I will focus on a few examples where the newly established transport laws have enabled us to make somewhat definitive statements about the role of fingering in explaining (or not explaining, in some cases) observations.
Understanding the surface metallicity of stars undergoing accretion of high-µ material from infalling planets is one such example. It is thought that, through the combined effect of type I migration in a protostellar disk (which can bring terrestrial planets close to the central star, [72]) and tidal interactions (which causes orbital decay of sufficiently close-in planets even in the absence of disk), planets could regularly fall into their host star even after the disk has disappeared [73]. The event is not expected to have a noticeable effect on the surface metallicity if the star has a deep outer convective zone, but could be important for stars with sufficiently shallow outer convection zones (or none at all). Planetary infall has thus been proposed as a possible mechanism that could explain the observed planet-metallicity correlation [34,74]. Vauclair [37] (see also [39]) correctly argued that this scenario would not work if mixing induced by fingering convection drains the metals into the interior on a short timescale compared with the time since the last infall event (for a given star). Using the mixing coefficient proposed by Traxler et al. [62] (see equation 39), which is valid in the limit where the density ratio is not too small (which is the case in these systems), I confirmed Vauclair's idea and demonstrated that any evidence for an infall event would disappear on a timescale of around 100Myr (see Figure 12). Since this is relatively short compared with the typical age of planet-bearing stars, we are left to conclude that the planet-metallicity correlation must be of primordial origin [38].
Fingering convection similarly affects the surface metallicity of WDs undergoing accretion from a planetary debris disk, and this effect should be taken account if one wishes to use the observed metallicities to derive the debris accretion rates [40][41][42].
Meanwhile, RGB stars are a good example of objects where a better understanding of fingering convection has made it harder to explain observations. Detailed spectroscopy of metal-poor stars by Gratton et al. [75] has revealed that the surface abundances of lithium and of elements participating in the CNO cycle changes noticeably along the RGB. A first sudden change occurs as expected during the first dredge-up event (i.e. when the outer convection zone penetrates most deeply into the star), see [76] for early work on the topic. A second much more unexpected change occurs once the convection zone has retreated, around the socalled luminosity bump that corresponds to the time where the hydrogen-burning shell enters the region that was previously mixed in the first dredge up, see Figure 13.
Eggleton et al. [49] noted that 3 He burning ( 3 He + 3 He → 4 He + p + p) is the dominant reaction in the cooler outer edge of the hydrogen burning shell, and locally lowers the mean molecular weight slightly. This can cause an inversion of the mean molecular weight gradient, but only once the shell has moved into a region that was previously homogenized by the dredge-up (see Figure 14). In the 3D simulations of Eggleton et al. [49], this inversion caused the development of a Rayleigh-Taylor instability, which mixed material between the hydrogen burning shell and the convection zone above. Charbonnel & Zahn [50] however pointed out that this inverse µ-gradient would first become unstable to fingering convection, rather  but not if C f = O(10). At the time, this was viewed as evidence in favor of the Ulrich prescription for fingering convection [9], but we now know that such a high constant C f is not supported by the numerical experiments presented in this lecture (see Figure 7 and also [60,77]).  In that case, an inverse µ gradient forms that drives fingering convection between the hydrogen burning shell (HBS) and the convection zone (CZ). This is a very frustrating state of affairs, since in all other respects the Charbonnel & Zahn model [50] provides a simple and elegant explanation for the data: fingering convection is triggered at exactly the right time and in the right place to cause the extra mixing needed -but its efficiency is too low. With this in mind, we dedicated the last 5 years (since 2013) trying to determine if another process, in combination with fingering convection, could increase the induced turbulent transport. This was at the heart of our attempts to establish whether thermo-compositional staircases could form through fingering in stars, but it really seems that they cannot (see Section 3.4 and [20]), at least spontaneously. It was also what prompted us to study added physics, such as rotation [23], magnetic fields [24] and shear [25], all of which are present in RGB stars. Of all these processes, magnetic fields appear to be the most promising. Preliminary results for fingering convection in a medium permeated by a vertical background magnetic field of amplitude B 0 obtained by Harrington and Garaud [24] suggest that for sufficiently large magnetic field strengths, where µ 0 is the permeability of the vacuum, and all the other quantities were defined earlier. This coefficient can be two orders of magnitude larger than one appropriate for non-magnetic fingering for B 0 as low as a few hundred Gauss, which are plausibly present in these stars [24]. More work is needed, however, to establish whether this result holds for arbitrarily aligned magnetic fields.

Oscillatory double-diffusive convection and layered convection 4.1 Traditional models of mixing by ODDC
As discussed in Section 1, the existence of regions that are both Schwarzschild-unstable but Ledoux-stable predates the discovery of ODDC by several years, and ODDC did not really become widely known in the astrophysical literature until much later. As such, there are many different prescriptions for mixing in semiconvective regions that have very little to do with the process of ODDC. These prescriptions are usually quite simplistic, either assuming that the region is adiabatically stratified, or purely radiative, or some interpolation between the two regimes, with some associated mixing model for chemical species.
It was only later that the first attempts to create models of semiconvective mixing based on the physics of ODDC were put forward [10,12,13,31,78]. Owing to its simplicity, the model of Langer et al. [12] is probably the most commonly-used one in stellar evolution codes today. It can be derived using the same arguments as the ones I presented in Section 3.1 to recover the Ulrich [9] and Kippenhahn et al. [11] models for fingering convection. If we simply assume that semiconvection acts as a turbulent diffusivity with coefficient D C = D semi , then by dimensional arguments D semi ∝ λd 2 =λκ T As discussed in Section 2.2.2, for sufficiently large wherek h is close to one. As a result, we can estimate D semi as where C L is a constant factor that can only be determined by comparison with experimental data. This expression is very reminiscent of the traditional fingering coefficient D fing (see equation 32), with R 0 replaced by R −1 0 . Using astrophysical notations that might be more familiar to this audience, this becomes which recovers the prescription of Langer et al. [12]. The coefficient C L is usually argued to be of order unity. Note that heat transport in this model is assumed to be negligible, so the background temperature gradient is the radiative one. Less well know perhaps is the model of Stevenson [31], who argued that the saturation of the unstable oscillatory modes of ODDC arises from a parametric subharmonic instability, in which smaller scale modes rapidly grow once the parent mode has reached a certain amplitude. He argues, using analogies with the geophysical literature, that the mixing coefficient should take the form 7 D which, for R −1 0 1, only differs from the Langer et al. proposal in the exponent applied to the inverse density ratio R −1 0 . Stevenson and his colleagues [31,78] were also the first to clearly argue for the existence of a distinct regime where layered double-diffusive convection takes place (though Spiegel [8] hinted at its possibility) instead of the small-scale turbulence considered so far, and to discuss where in parameter space it could take place.
Indeed, as discussed in Section 2.2.2, the range of linear instability for ODDC in geophysics (and in laboratory experiments) is very small, so ODDC almost never derives from it. Instead, it is excited from a subcritical branch of instability. As a result, none of the existing laboratory experiments on ODDC show the presence of unstable gravity waves, but instead, all take the form of layered double-diffusive convection (layered convection, for simplicity) [19,[79][80][81]. Layered convection is also ubiquitously found in nature on Earth in the combined presence of unstable temperature gradients and stable salt gradients, such as in volcanic lakes [82], in the arctic ocean [83], and under the ice shelf in the antarctic [84]. The temperature and salinity profiles associated with layered convection take the form of a thermohaline staircase similar to fingering staircases discussed in Section 3.4, although the overall gradients now have the opposite signs.
An entirely different class of astrophysical semiconvection models therefore exists that assumes the presence of layered convection [13,78,[85][86][87]. Despite their differences, these models generally start from the same assumption, namely that the thermocompositional staircase is in equilibrium. This implies that the fluxes of T and C through a layer must be equal to the fluxes of T and C through the adjacent interfaces on either side, otherwise the interfaces and/or the layers would have to evolve on a rapid timescale 8 . As a result of this assumption, one can easily build a transport model from two ingredients only: the first is a prescription for the temperature flux within a convective layer (which is not negligible in layered convection), and the second is a prescription for the ratio of temperature and compositional fluxes across an interface. Once these two quantities are known, the fluxes of heat and composition through the entire staircase can easily be computed.
The turbulent temperature flux through individual layers can be put forward on dimensional grounds, to be where Nu T is a potential temperature Nusselt number, defined in Table 1; going from the second to the third expression above uses that definition. Traditional geophysical models of overturning convection between solid boundaries separated by a distance H usually argue [64,[88][89][90][91] that for very turbulent flows the Nusselt number should be equal to some power of the so-called Rayleigh number Ra = |N 2 T |H 4 /κ T ν, which is the ratio of the buoyancy force driving convection, to the viscous force that damps fluid motions. In layered convection in stars, however, the viscous force is thought to be negligible, and the relevant Rayleigh number is instead, where H L is the layer height. As such, models of layered convection in stars traditionally have where the pre-factor C semi , and the power a, vary between models (for instance a = 1/4 in [13] and a = 1/3 in [85,86]). We see that F T depends rather sensitively on the layer height, which must also be specified by the model. Without further justification, H L is usually taken to be some fraction of the pressure scaleheight, to be specified by the user. Note that unless H L is very small, Nu T is much larger than one, so the diffusive contribution to the heat transport is in this case negligible. The second ingredient of these models is a parametrization for the non-dimensional ratio of the temperature flux to the compositional flux, γ −1 tot (see Table 1). In this case, all traditional astrophysical models of layered convection [13,78,85,86] agree with one another and are based on laboratory experiments [79,92] and theory [80] of layered convection in salt water. In this theory, the interface is assumed to be entirely diffusive, so the ratio of the interfacial fluxes is equal to the ratio of the diffusive fluxes, which in turn depends on the respective temperature and salinity gradients across the interface. With further assumptions, Linden & Shirtcliffe [80] arrive at the conclusion that (with no proportionality constant between the two expression). This prediction is consistent with many of results obtained in laboratory experiments, and is generally considered to be correct in the geophysical literature.
This flux ratio can then be used to compute the compositional flux across the staircase given the temperature flux, as in which in turn implies that the compositional mixing coefficient for layered semiconvection can be written as where the Nusselt number is given in (83) and γ −1 tot = τ 1/2 . Again, unless the layer height is very small, this quantity is usually much larger than the microscopic diffusion coefficient κ C , which can be neglected.

Numerical simulations of ODDC
As in the fingering case, none of these models had ever been tested in the low Prandtl number regime more appropriate of stellar astrophysics until recently. In particular, whether ODDC takes the form of small-scale wave-like turbulence or layered convection remained unknown. Early simulations of ODDC/semiconvection in a two-dimensional vertically bounded domain were presented by Merryfield [57] (see also [93]), but the results were limited in scope by the resolution affordable at the time. The first three-dimensional direct numerical simulations of ODDC at low Prandtl number were presented by Rosenblum et al. [94], using the PADDI code [63] with the model setup described in Section 2.1, and solving equations (11)-(13) with the − sign in the temperature and composition equations. In this paper, we demonstrated that both layered and non-layered outcomes are possible, depending on the local thermocompositional stratification (as measured by the inverse density ratio R −1 0 ). This was later confirmed by Mirouh et al. [22], who performed a more comprehensive exploration of parameter space. These two possible outcomes are illustrated in Figure 15.
For high inverse density ratios, i.e. for systems that are more strongly stratified, we have found that ODDC is excited as expected from the linear theory described in Section 2.2.2, and saturates into a state of weak wave turbulence. The fastest-growing modes are, as expected, elevator modes. However, once they reach a certain amplitude, nonlinear interactions causes a transfer of energy to modes that have a higher vertical wavenumber (see Figure 15b). This is at least qualitatively consistent with the idea put forward by Stevenson [31]. The turbulent fluxes in that regime are fairly weak, and decrease with increasing R −1 0 , and decreasing Pr and τ (see below).
For low inverse density ratios, on the other hand, the initial state of wave-like turbulence always transitions to a layered state. The initial height of the convective layers is fairly small, of the order of a few tens of d, but the layers then always merge until a single one is left in the computational domain (the mergers cannot proceed beyond that owing to the periodicity of the boundary conditions). The initial formation of the layers, and each subsequent merger, is accompanied by a substantial increase in the turbulent fluxes, suggesting that the latter indeed depend on the layer height.
In view of these results, it is clear that any model of ODDC/semiconvection in stars should include a way to determine which regime is expected (layered or non-layered), a prescription for transport in the layered regime and a prescription for transport in the non-layered regime [22].

Transport in non-layered ODDC
Starting with the latter, we can use the simulations of Mirouh et al. [22] to extract turbulent fluxes in non-layered convection 9 , and compare them to the predictions of Stevenson [31] and Langer et al. [12] (which are for the non-layered regime). The turbulent diffusion coefficient D semi is extracted from the simulations as usual by computing − wC /C 0z = R 0 κ T ŵĈ . As a side note, it is easy to show using arguments very similar to those put forward in the fingering regime (see Section 3.2), that both ŵĈ and ŵT have to be positive in a statistically stationary state. As such, D semi is indeed positive.
The results are shown in Figure 16a. In the model of Stevenson [31], D semi /κ T ∝ R 2 0 while in Langer et al. [12] D semi /κ T ∝ (R −1 0 − 1) −1 . The data clearly favors a model where D semi /κ T is a power law function of R −1 0 − 1 and the slope is roughly equal to −1.5, which is somewhere in between these existing models, at least for large enough R −1 0 . We also see that these models fail to account for the stabilization of the system (with D semi → 0) that occurs when R −1 0 approaches the marginal stability threshold (Pr + 1)/(Pr + τ). In an attempt to create a new model for mixing by ODDC in the non-layered regime, we (i.e. me and my students) have tried various approaches, including an extension of the Brown et al. model to the ODDC case, and various attempts at weakly nonlinear theory to capture the nonlinear evolution of the unstable gravity waves. None of these approaches were able to model the variation of D semi with all input parameters (Pr, τ, and R −1 0 ) satisfactorily. This does not mean that such a theory cannot be developed -only that we have been unable to do so to date. For lack of a more theoretically-grounded model, Mirouh et al. [22] put forward an empirical parametrization of the existing data that seems to capture adequately the variation of both compositional and temperature fluxes across parameter space. In this parametrization, the non-dimensional turbulent temperature fluxF T = ŵT is proposed to beF where R −1 oc = (Pr + 1)/(Pr + τ) is the marginal stability threshold for ODDC, and is the reduced inverse density ratio associated with ODDC. This parameter, as in the fingering case, maps the instability range R −1 0 ∈ [1, R −1 oc ] to the interval r o ∈ [0, 1]. A comparison of (68) with data from [22] is shown in Figure 16b, suggesting that this models is adequate. We also note thatF T rapidly decreases with Pr and τ, to the extent that for any reasonable stellar valueF T is either smaller than or negligible compared to the diffusive flux (which is one in these units). We therefore confirm that the turbulent heat flux is negligible in comparison with the diffusive flux for non-layered ODDC in stars.
To model the compositional flux, we could use a similar empirical parametrization, starting for instance from the scaling with R −1 0 − 1 discovered in Figure 16a; this may be an interesting idea to pursue in the future. Generally speaking, the data suggests that D semi rapidly decreases with increasing stratification. Moll et al. [95] in fact suggested that compositional mixing in the non-layered regime (which takes place at large R −1 0 ) is negligible compared with pure diffusive transport at stellar parameters. This remains to be confirmed.
As an alternative approach, in Appendix A.3. of Mirouh et al. [22] we extended the work of Schmitt [14,96] (also used by Brown et al. [21]), who noted that it is possible to estimate the ratio of the turbulent temperature flux to the turbulent compositional flux, named γ turb (see Table 1) in double-diffusive systems simply from linear theory, even though nonlinear arguments are needed to estimate each of them individually. Indeed, using non-dimensional fluxes, we can write γ turb = ŵT / ŵĈ . Assuming thatŵ,T andĈ are dominated by the fastest-growing linearly unstable modes, we can replace them by their proposed ansatz from linear theory (see Section 2), to get In this expressionλ andk h are the growth rate and wavenumber of the fastest-growing elevator modes for the desired set of parameters Pr, τ and R −1 0 , computed from linear theory (see Section 2). From the linearized versions of equations (11)-(13), we have the following system of equations relating the amplitudesw,T andC of each mode.
where as usual + refers to the fingering regime and − refers to ODDC. Solving forT andC, we getT If we finally assume, without loss of generality, that the amplitudew is real, then we combine equations (70) and (72) to get the following prediction for the turbulent flux ratio γ turb [22]: where we have writtenλ =λ T + iλ I . Conveniently, the ± signs cancel out, so this expression can be used for both fingering convection and ODDC. Even more conveniently, the unknown amplitude of the modew cancels out entirely from that expression, so γ turb does not contain any unknown parameter. A comparison of this theory with the simulation data (in terms of γ tot , rather than γ turb ) is presented in the next section, and shows reasonably good (but not perfect) agreement between the two. This expression can finally be used to compute the dimensional mixing coefficient with ŵT given by (68). One may rightfully argue that this method is much more complicated than proposing a simple empirical formula for D semi based on the data shown in Figure 16, which would be more practical for stellar evolution purposes. I do not disagree, especially since we find that in most instances D semi κ C at stellar parameters, making the additional effort of computing γ turb (which involves solving a cubic equation) somewhat pointless. It does however have the advantage (in my opinion) of being based on solid theoretical arguments for why the turbulent flux ratio γ turb should be of a particular form, and can in addition be used to establish a criterion for layer formation, as shown below.

Criterion for layer formation
As discovered by Rosenblum et al. [94] and Mirouh et al. [22], simulations of ODDC at a low inverse density ratio always spontaneously evolve into a state of layered convection, suggesting that the mechanism leading to layer formation must be very robust. In Section 3.4, we learned about the γ-instability as a mechanism for layer formation in oceanographic fingering convection, and found that it should occur whenever the total flux ratio γ tot is a decreasing function of the density ratio R 0 . It is therefore natural to wonder whether a similar mechanism may be at play here.
As it turns out, the mechanism is not just similar, it is indeed exactly the same. To see this, simply note that the horizontally-averaged temperature and composition equations which were used as the starting point for the development of the γ−instability theory in Section 3.4 (see equation 45) are exactly the same in fingering convection and in ODDC, since the horizontal average of the vertical advection terms containing ±w are identically zero (for mass conservation). As a result, Rosenblum et al. [94] demonstrated that the criterion for layer formation in ODDC is the same as in fingering convection: layers are expected to form when γ tot is a decreasing function of R 0 , or equivalently, when γ −1 tot is a decreasing function of R −1 0 . This theory can be tested from the extended dataset of Mirouh et al. [22], and some of the results are presented in Figure 17a. All symbols in this figure show the measured values of γ −1 tot extracted from the data, by computing in the non-layered phase (this expression is the non-dimensional version of the one given in Table 1). Large symbols in this figure denote simulations that eventually transitioned into layered convection. We see that this indeed happens whenever γ −1 tot is a decreasing function of R −1 0 , confirming the role of the γ-instability as suspected.
A more quantitative way of testing the theory is to compare the predicted and observed growth rates of the layering modes in simulations. Let's take for instance a simulation presented by [94], with Pr = τ = 0.3 and R −1 0 = 1.2. In this simulation, 4 layers were first observed to form, though by inspection of the spectral power in the horizontally-averaged density field (see Figure 17b) we see that the 3-layer mode is also growing. In order to compute the theoretical growth rateΛ of the layering modes, we first need to estimate the parameters of the quadratic dispersion relation (54). We can easily measure the potential temperature Nusselt number Nu 0 and the total flux ratio γ 0 in the non-layered phase preceeding layer formation in the simulation. We can also, by running simulations at slightly different R −1 0 , compute A Nu = R 0 for Pr = τ = 0.3 and R −1 0 = 1.2. Using these parameters and coefficients in the quadratic equation (54), we then predict thatΛ (K) 0.47K 2 (77) for this simulation. We can compare the temporal evolution of the spectral power in density for each growing mode with this prediction, knowing their wavenumber isK = k n , where k n is the wavenumber with n layers in the domain. We see in Figure 17b that the growth rate of the 4-layer mode is over-predicted by the model, but that of the 3-layer mode is adequately captured. As discussed in Section 3.4, the fact that the mean-field theory fails when the mode wavenumber is large is not surprising given the ultraviolet catastrophe inherent to the mean field model, but it is very encouraging to see that it works well for the 3-layer mode. Similar comparisons were later made by Mirouh et al. [22], with equally satisfactory success. This completes the validation of the γ−instability theory for layer formation in ODDC.
If one believes the Mirouh et al. [22] prescriptions for the turbulent heat flux and turbulent flux ratio described above, it is easy to predict when layers are expected to form in stellar ODDC. Indeed, from (75) and the definition of γ turb = ŵT / ŵĈ we have  Figure 17. Left: The symbols show the total inverse flux ratio γ −1 tot measured from simulation data in non-layered ODDC for various Pr and τ (see [22] for a more extended dataset). Larger symbols denote simulations that ultimately become layered. The solid lines are predictions for γ −1 tot from the model of Mirouh et al. [22]. Figure adapted from [1]. Right: The solid lines show the growth of the spectral power in the horizontally-average density field for various layering modes with wavenumber k n (where n is the number of layers in the domain) for a simulation with Pr = τ = 0.3 and R −1 0 = 1.2 [94]. The 3-and 4-layered modes grow until they reach an amplitude that causes inversions in the density profile (black line), at which point layered convection begins. The 1-and 2-layered modes do not initially grow but emerge out of subsequent layer mergers. The dotted lines show the corresponding prediction from the γ−instability theory for the growth of the 3-layered and 4-layered mode.
We can then use (68) and (73) to predict γ −1 tot as a function of R −1 0 for any parameter set (Pr, τ). A comparison of this prediction with existing data is shown in Figure 17a, and reveals very good agreement (see [22] for a comparison with a more extended dataset), giving us some confidence in the model. Using this method, we can determine when γ −1 tot is a decreasing function of R −1 0 for any Pr and τ in the stellar parameter regime. The results were first presented by Mirouh et al. [22] and are shown in Figure 18, which is a contour plot of the critical value of R −1 0 below which layering can occur, as a function of both Pr and τ. Generally speaking, we see that when Pr > τ 1/2 the layering threshold is more-or-less independent of τ and is roughly equal to Pr −1 . This is not a particularly relevant situation in stellar astrophysics, however. On the other hand, when Pr < τ 1/2 , then the layering threshold is roughly equal to τ −1/2 (with a weak dependence on Pr). For stellar interiors with Pr and τ of order 10 −6 , we therefore expect layers to spontaneously form whenever the inverse density ratio R −1 0 < O(10 3 ), implying that all cases of semiconvective zones near a convective core should lie in the layered regime. For planetary interiors or degenerate regions of stellar interiors with Pr and τ of order 10 −2 , we expect layers to spontaneously form whenever the inverse density ratio R −1 0 < O(10) [22].

Mixing in layered convection
Having established when layered convection is naturally expected, we now turn to the question of modeling turbulent transport induced by the process. Figure 15a revealed that the total heat and compositional flux through layered convection depends on the layer height, and increases every time a merger occurs. To study this more quantitatively, Wood et al. [97] (see also earlier preliminary work by [94]) re-analyzed the simulations of Mirouh et al. [22] that became layered, and also presented new simulations that were produced using the same model setup and code, but in much larger domains. They measured the average turbulent temperature and compositional fluxes in the domain as a function the average layer height. With this information, the Rayleigh number Ra , the potential temperature Nusselt number Nu T and the total flux ratio γ −1 tot can be computed and compared with the models introduced in equation (63) for layered convection.
The results for existing simulations with varying Pr, τ and R −1 0 are shown in Figure 19. Figure 19a reveals that Nu T is indeed a power law function of Ra , with: with C semi ∼ O(0.1) when fitted to the existing data. The exponent of that power law therefore rules out [13] but is consistent with [85,86]. The pre-factor C semi seems to have an additional weak dependence on both R −1 0 and τ, which could not be determined without much more demanding numerical simulations [97]. Indeed, we just saw that layered convection only spontaneously emerges when R −1 0 is smaller than the layering threshold (which is ∼ τ −1/2 ). As such, running simulations of layered convection for a wider range of R −1 0 requires using much smaller values of τ, which has so far been computationally prohibitive. Figure 19b tests the second assumption of traditional models [13,78,85,86], namely that γ −1 tot = τ 1/2 [80]. We see that this prescription does not hold: γ −1 tot is generally larger that τ 1/2 (see also the results of Moll et al. [33]), having instead a complicated dependence on Pr, τ and R −1 0 , and perhaps even on the layer height (though there are hints that γ −1 tot may become independent of H L for sufficiently tall layers). The discrepancy between these numerical results and the Linden and Shirtcliffe [80] model for the interfacial flux ratio is easy to understand. The interfaces in their model (and more generally in the high Prandtl number laboratory experiments of layered convection for which the model was created) are quiescent so the interfacial transport is purely diffusive. However, the interfaces in our low Prandtl number simulations are quite turbulent (see Figure 15). This example nicely illustrates the risk of applying geophysically-derived parametrizations to stellar astrophysics without verifying them first with low Prandtl number numerical experiments.
Based on their limited data, Wood et al. [97] proposed the following empirical parametrization for the compositional mixing coefficient: Note, however, that combining this with (79) would imply that (neglecting very weak powers of Pr and Ra and assuming that Nu C , Nu T 1). This is not really consistent with Figure 19b, where it appears that γ −1 tot may be decreasing instead of increasing with R −1 0 . This discrepancy most probably comes from the fact that neither (79) nor (80) correctly account for the weak dependence of Nu T and Nu C (and in turn D semi ) on   R −1 0 . This problem can only be resolved by running simulations for a wider range of R −1 0 , which is computationally very challenging as discussed above.
Ultimately, however, the main uncertainty in modeling both heat and compositional fluxes in layered convection is their dependence on the layer height (i.e. with F T,C ∼ H 4/3 L ). In all the simulations performed to date, the layers were ultimately seen to merge until a single one remains in the domain. This poses two problems. First, it has prevented us from proposing a model for what the layer height may be in a star, since it is not clear what mechanism (if any) eventually stops the merger process. Secondly, these results seem to challenge the original assumption of an equilibrium staircase, suggesting instead that layered convection may be a short temporary phase in the evolution of a star, whose duration depends on the layer merger rate.

Conclusions for now
Our numerical results and theoretical investigations have confirmed that there are two regimes of ODDC / semiconvection: one that is layered and one that is not layered. The layered form arises spontaneously from the γ−instability [22,69,94] at sufficiently low inverse density ratio. The threshold for the layering instability to occur can be computed in the manner described in Section 4.4 [22], and is approximately given by R −1 0 τ −1/2 at stellar (and planetary) parameters. Note that layered convection can also exist for R −1 0 beyond that threshold, but only if layers are seeded, i.e. if they are present in the initial conditions applied to the problem (see [57], [93], [33]). Whether this regime is relevant for stars or not remains to be determined, since one would have to specify why these layers should be present in the first place.
Both turbulent heat and compositional transport in layered convection are very efficient, and can be computed as follows [97]: with where Ra depends on the layer height as H 4 L (see equation 63), and C semi is of the order of 0.1, with some possible weak dependence on τ and R −1 0 that remains to be determined. What the layer height would be in a star, or whether the concept of an equilibrium staircase even applies in the first place, remains to be determined. Meanwhile, the turbulent mixing where γ −1 tot is of order one (but smaller than one), with some possible weak dependence on Pr, τ and R −1 0 that also remains to be determined. In the non-layered regime, by contrast, we have demonstrated that the turbulent heat flux is negligible, and the turbulent compositional flux also appears to be negligible, at least for non-degenerate regions of stellar interiors where Pr and τ are asymptotically small. At more moderate Pr and τ, the turbulent compositional flux could be significant, and can be determined from the model of Mirouh et al. [22] presented in Section 4.3, or by fitting their data empirically (which remains to be done).

Applications to stellar astrophysics
The impact of these newly developed theories of mixing by ODDC (layered or non-layered) on stellar evolution have not really been studied in detail yet, so this section will be very limited in scope.
As discussed in Section 2.4, semiconvective regions are most commonly found just outside the convective cores of MS stars, in two mass ranges: intermediate-mass stars with M between about 1M and 3M , and higher-mass stars with M > 10M . In both cases, the typical inverse density ratios are well within the layered regime, with R −1 0 τ −1/2 . The turbulent transport of heat and composition can be substantial and depends on the unknown heights of the layers (which need not be constant in space nor time). With all this uncertainty, the best one can do (for now) is to fix the layer height to some fraction of a pressure scaleheight, and qualitatively ascertain what the effects of layered convection on the star's evolution are as a function of the assumed layer height .
It has long been known that one of the main effects of semiconvection in these stars is to prolong their lifetime on the MS by increasing the amount of hydrogen available for fusion reactions. Consequently, it also increases the size of the helium core at the end of the MS. Moore and Garaud [98] attempted to quantify the effect, to determine more specifically how the helium core size of intermediate-mass stars varies with the assumed semiconvective layer height. To do so, Moore implemented both heat and compositional transport prescription by layered ODDC in MESA. The compositional transport is easy to implement because it merely requires computing the mixing coefficient D semi with equation (67). Including a model for the turbulent heat transport on the other hand is much more complicated, since MESA computes the temperature profile within the star given the heat flux, instead of computing the heat flux given the temperature gradient (which is the way in which equation 82 is written). As such, including heat transport by layered semiconvection involves the numerical solution of a fourth-order polynomial in the quantity ∇. Details of the implementation are given in Section 3.1.2 of [98].
The conclusions of our investigation on intermediate-mass stars are remarkably simple. We found that unless the height of the layers was assumed to be unphysically small, compositional mixing by layered convection was so efficient that it would erase any existing compositional gradient on a short timescale compared with the stellar evolution timescale. As a result, the semiconvective region rapidly becomes fully convective, and the star's core grows to be the size one would have predicted much more simply using the Schwarzschild criterion and ignoring semiconvection altogether. This is illustrated in Figure 20, which shows the predicted (fully) convective core size of a 1.3M star, under four different scenarios: ignoring semiconvection and with either the Schwarzschild or Ledoux criteria to compute the core sizes, or, including layered convection and assuming that the layer height H L is either 10 −6 or 10 −14 times the pressure scaleheight H p . We see that with H L = 10 −6 H p the model already behaves almost as if the Schwarzschild criterion had been used instead, and it is only for vanishingly small layer heights (H L = 10 −14 H p ) that a substantial semiconvective region survives above the core throughout the star's evolution. These results are consistent with recent asteroseismic determinations of the convective core sizes in MS stars [99,100], which reveal them to be even larger than what stellar evolution predicts using the Schwarzschild criterion. This implies not only that semiconvection is irrelevant (as we predicted), but also that substantial overshoot must be present. We therefore conclude that ignoring semiconvection and using the Schwarzschild criterion in these stars is a good approximation for stellar evolution purposes (for intermediate mass stars), which is particularly convenient given how technically difficult the implementation of layered convection is (for heat transport) in standard stellar evolution codes.
Finally, I would like to add that MESA, to my knowledge, is the only stellar evolution code to date that includes the complete Wood et al. [97] prescription for transport by layered convection for both heat and composition in semiconvective zones. I encourage readers to use it to see what other effects layered convection may have, in particular for more massive stars (M > 10M ), or for stars in more evolved stages of evolution.

What next?
In the previous Sections, I have attempted to summarize what is known so far about doublediffusive instabilities and double-diffusive convection at low Prandtl number, which is the relevant limit for astrophysical fluids. However, there are a number of limitations of the model we have used to arrive at these results, which I would be remiss not to mention. Much also remains to be done in terms of modeling double-diffusive convection more realistically, including for instance the effects of rotation, magnetic fields, shear, and other relevant processes. I now discuss both of these in turn.

Model limitations
The numerical model used in the simulations presented in this lecture uses the Spiegel-Veronis-Boussinesq [17] approximation, and triply-periodic boundary conditions, so it is important to determine when the limitations of this setup may begin to invalidate the results obtained.
The Spiegel-Veronis-Boussinesq approximation is valid a long as two conditions are satisfied: (1) that the typical fluid velocities are much smaller than the sound speed c s and (2) that the region modeled is much smaller than a pressure scaleheight H p . In Section 2.2, we learned that the typical size of basic double-diffusive structures (prior to the formation of layers, if they form) is of order 10d = 10(κ T ν/|N 2 T |) 1/4 and provided some estimates for that quantity. We also saw that the structures remain small even in the nonlinear regime (again, prior to layer formation), so 10d is a good estimate for the typical eddy size in small-scale homogeneous double-diffusive convection (of both fingering and ODDC type). In terms of the typical velocity of fingering convection and ODDC, we saw that it can be estimated using λd, where the growth rate λ is at most equal to N T , but usually much smaller (depending on the regime considered). It is clear that well within the interior of a star, we always have 10d H p and N T d c s . Closer to the stellar surface, however, especially for intermediate and high-mass stars, κ T increases rapidly, and so will d. Meanwhile, the pressure scaleheight H p and the sound speed c s both decrease, so it is conceivable that fingering convection or ODDC in regions sufficiently close to the surface of these stars may need to be modeled more realistically including the effects of compressibility.
Once layers form in ODDC, the typical size of turbulent structures increases to be the entire layer height, and the turbulent velocities also increase to take typical convective (rather than double-diffusive) values. If the layers grow to a substantial fraction of a pressure scaleheight, the effects of compressibility should therefore also be taken into account.
Finally, note that the Boussinesq approximation also assumes that the diffusivities are constant within the domain. This could be a limitation of the model when ν, κ T or κ C depend sensitively on temperature or composition.
Another set of caveats of the results comes more specifically from the use of the triplyperiodic boundary conditions for the perturbations in our DNSs. The first is relatively well understood, and pertains to the relative size of the double-diffusive structures modeled compared with the computational domain size. If the structures are much smaller than the domain size (in all directions), then the latter has no effect on the double-diffusive dynamics. This is usually the case in the homogeneous regime of fingering convection and ODDC. On the other hand when the double-diffusive structures approach the domain size, it is quite likely that the results become dependent on the boundary conditions. This is the case for instance in layered convection, when there is only a single layer left. This is also the case of rapidly rotating double-diffusive convection and magnetized fingering in the very strong field regime, where the dynamics become invariant along the rotation axis and the magnetic field, respectively. When this happens, the measured fluxes can be artificially enhanced [23,24,29]. To control this problem, it is important to run additional DNSs in increasingly large domains until a point where the results become independent of the domain size. In all the simulations presented in these lectures, this was always verified.
A second problem associated with triply-periodic boundary conditions is more subtle, and requires care in interpreting the turbulent flux data. In our computational model setup, the background potential temperature and composition gradients are fixed, and the developing instabilities drive turbulent heat and composition fluxes. Crucially, there is nothing in the system that limits these fluxes beyond what drives them. However, in a real star the situation is somewhat different because the total heat flux that needs to be transported through a given shell of material is fixed by the nuclear reaction rates taking place within or below. In that case, we anticipate that (on long timescales) the double-diffusive dynamics may adjust themselves to carry the required flux, rather than the flux being solely controlled by the double-diffusive dynamics. The problem is not significant when the double-diffusive dynamics transport a negligible amount of heat (which is the case of homogenous fingering convection and ODDC), but needs to be born in mind for layered convection and magnetized fingering convection with strong fields, where heat transport can be substantial (see below). These types of problems are then perhaps best studied in a model setup that can impose a fixed heat flux through the domain.

More realism
In the past few years, our group has begun to investigate the effects of other physical processes on the dynamics of double-diffusive convection, including those of lateral gradients of potential temperature and composition [26], rotation [23,29] magnetic fields [24], and shear [25]. Given the dramatic increase in the dimensionality of parameter space that needs to be explored, our results are in most cases preliminary, and much remains to be done. This section briefly summarizes the state of the field as of May 2019.

Horizontal gradients
The effects of horizontal gradients have been investigated in the context of fingering-like stratifications only [26] (the ODDC-like case has never been looked at, to my knowledge), and affect the results described in Section 3 only in the limit of very strong stratification, i.e. for density ratios approaching marginal stability or exceeding it. Horizontal gradients greatly increase the region of parameter space that is unstable [5,26], to the extent that almost any inverse µ gradient can cause vertical mixing. However that instability is slowly growing, and does not appear to cause very strong enhancements in the turbulent mixing coefficient above the purely diffusive value. This conclusion was obtained from simulations at very moderate Prandtl number and diffusivity ratio, and it remains to be seen whether it holds for more stellar-like parameter values.

Rotation
The effects of rotation have been investigated in preliminary work by Moll and Garaud [29] for ODDC, and Sengupta and Garaud [23] for fingering convection. In both cases, rotation does not affect the overall instability range, nor the growth rate of the fastest-growing modes of instability. Whether rotation has a weak or strong effect on the nonlinear saturation of double-diffusive instabilities can be determined by computing the turbulent Rossby number, Ro = u rms /Ωd where u rms is the typical velocity of turbulent eddies, and Ω is the rotation rate. Rotation is important for Rossby numbers of order one or much smaller than one, but is negligible when Ro 1. At moderate Pr and τ (appropriate for giant planets, or degenerate regions of stellar interiors), the turbulent velocities can be estimated roughly using κ T /d, so the Rossby number is Ro = κ T /Ωd 2 = Pr −1/2 |N T |/Ω [29]. With that assumption and the fact that Pr ∼ 0.01 in these objects, we see that the effects of rotation will only be negligible if Ω 10|N T |, roughly speaking. In stellar interiors, the turbulent velocities associated with the basic instability are typically much smaller, of the order of (κ T /d)(Pr/R 0 ) 1/2 for fingering convection and (κ T /d)(PrR 0 ) 1/2 for non-layered ODDC. As a result, the Rossby number is of the order of Ro = R −1/2 0 |N T |/Ω for fingering convection and Ro = (R −1 0 ) −1/2 |N T |/Ω for nonlayered ODDC. In other words, at very low Pr and τ and for sufficiently strongly stratified systems (i.e. large R 0 for fingering, and large R −1 0 for ODDC), the Rossby number can be of order one or even smaller and rotation could be important. For layered ODDC, the situation is a little different since the estimate for the turbulent velocities u rms must be replaced by an estimate for convective velocities, that depends on the assumed layer height.
Most of our numerical simulations of rotating double-diffusive convection to date have been performed at the poles (where the rotation axis is aligned with gravity), with a few cases only looking at mid-latitudes. We have found that weak rotation (with Ro 1) tends to reduce the vertical compositional flux slightly compared with non-rotating simulations, in both fingering and ODDC. By contrast, when the Rossby number is of order unity, we found that rotation causes the formation of large-scale vortices that span the entire domain in both vertical and horizontal directions 10 . These vortices have a tendency to cause a segregation of the compositional field in their core, either expelling it or concentrating it (see Figure 21) depending on the initial conditions selected. Denser cores drive stronger downward velocities, while lighter cores drive stronger upward velocities. This increase in the correlation between w and C causes a substantial increase in the vertical compositional flux compared with a non-rotating simulation at otherwise similar parameters. The effect is very substantial, and could potentially be invoked as a solution to the missing mixing problem in RGB stars (see Section 3.6) [23]. However, why these vortices form, and whether they would also form in stars (rather than in idealized simulations), remains to be determined. In particular, based on our results, and similar results obtained for rotating Rayleigh-Bénard convection [101,102], it is not clear that they would form at lower latitudes [29], or in computational domains that have different horizontal dimensions (i.e. L x L y ) [103]. . This is the compositional field perturbation, which shows that the vortex core is low C, and flows upwards. Figure adapted from [23]. Right: Illustrative snapshot of a large-scale vortex formed in an ODDC simulation at Pr = τ = 0.1, and R −1 0 = 1.25, with Ta * = 10. This time the vortex core is high C, and flows downwards, but has a more complex radial structure that remains to be explained. Figure adapted from [29], Finally, for rotationally dominated systems (with Ro 1), large-scale vortices do not form. Instead, the dynamics of both fingering and ODDC are dominated by narrow tubular structures that are invariant along the axis of rotation. In that case, the periodicity of the domain becomes a limitation, and it is not clear that the numerical results obtained are reliable 10 Note that they have only been seen in polar cases so far (see previous section). More remains to be done, using perhaps a different formalism (or much taller periodic domains) to characterize ODDC and fingering convection in that limit.

Magnetic fields
As discussed earlier, Charbonnel and Zahn [104] and Stevenson [31] investigated the linear stability of magnetized fingering convection and magnetized ODDC, respectively, but very little is known so far about the effects of magnetic fields on double-diffusive convection in the nonlinear regime.
Beyond linear theory, the first numerical study of magnetized fingering convection has just been completed by Harrington and Garaud [24], in a model that assumes the presence of a vertical 11 large-scale background field of amplitude B 0 . The impact of the field can be characterized by a new non-dimensional number, namely which is the square of the ratio of the Alfvén velocity associated with that field, to the selected unit velocity κ T /d (see Section 2.2). Note however that since κ T /d is not always a good estimate for the actual fingering velocities (see Section 3), especially at very low Prandtl numbers, H B 1 does not necessarily imply that the magnetic field is negligible (see below for more on this topic). We have run a number of DNSs of fingering convection that include a vertical magnetic field, varying the quantity H B , with otherwise fixed parameter values. We have found that vertical transport can be increased by orders of magnitude compared with non-magnetic simulations, because the field stabilizes the fingers partially against the parasitic shear instabilities that normally destroy them (see Section 3). As a result, larger vertical velocities can be achieved before the shear instability causes saturation (see Figure  22), which in turn enhances the vertical fluxes. Our results can be quantitatively explained by a simple extension of the Brown et al. [21] model, that contains no additional free parameter. Recall that the mixing coefficient D fing can be estimated using equation (42), where K B 1.24 and whereŵ f is the anticipated vertical velocity within the finger at saturation obtained by matching the finger growth rateλ to the parasitic shear instability growth rateσ. The latter decreases in the presence of a magnetic field [105], because the field rigidifies the fingers. This effect can be captured using the approximate expression (see [24] for detail): This equation recovers the hydrodynamic limit studied by Brown et al. [21] when H B = 0, as required. It also shows that when H B 0 the fingering elevator modes are unaffected by parasitic shear instabilities until their vertical velocities reach a threshold value ofŵ f = √ 2H B (which corresponds to energy equipartition). Beyond that threshold, the finger velocities continue to grow untilσ C Hλ , where C H 1.66 is not a free parameter, but instead, is fixed by the requirement that our theory recovers the hydrodynamic limit when H B = 0 (see [24] for detail).
All that remains is to solve the algebraic equation C Hλ = 0.42ŵ fkh (0.5 − H Bŵ −2 f ) 2/3 for w f numerically, then substitute the result in equation (42) to compute the mixing coefficient D fing . Note that for large H B , we find thatŵ f √ 2H B so using this expression in (42) yields which recovers the aforementioned equation (57) when K B 1.24, τ 1, and if we use the approximationλ √ Pr/(R 0 − 1) (see equation 18). Figure 24 compares the new model for magnetized fingering convection with results of our DNSs including the vertical field, at moderate parameter values (Pr = τ = 0.1, R 0 = 1.45), and confirms that it correctly predicts the variation of D fing with magnetic field strength as H B increases. Applying the model to more stellar-like conditions (where the fingering growth rate is much smaller), we predict that magnetic fields become important for H B ∼ 10 −4 in WDs, and 10 −8 in RGB stars. With a very moderate magnetic field of the order of 300G, D fing increases by two orders of magnitude compared with the non-magnetic case for both types of stars (see Figure 24). This strongly suggests that magnetic fields should be taken into account when modeling fingering convection in stars, and could be the answer the RGB star abundance conundrum described in Section 3.6.
D fing κ C Figure 23. D fing /κ C in simulations (symbols) and models (lines). The simulations were run at Pr = τ = 0.1 and R 0 = 1.45, and the model corresponding to these parameters is the black line. The red and blue lines are predictions for RGB and WD stars respectively, see [24] for detail. The shaded areas are values of H B corresponding to B 0 ∼ 100G. Figure adapted from [24].
Of course, this conclusion remains to be tested in the more general case where the magnetic field is inclined with respect to gravity. Preliminary results suggest that transport is also increased in that case, but new effects also appear (such as a dynamo, see [24] and Harrington's MS thesis), whose impact needs to be better understood.

Shear
The final physical process whose impact on stellar double-diffusive convection had not yet been investigated yet is shear, though numerical and experimental studies of sheared doublediffusive convection in geophysics exist, such as for instance [106][107][108][109] for the fingering case, and [30,110] for ODDC.
Although sheared ODDC in stars has not been studied yet, the effect of a moderate shear on fingering convection is now relatively well understood [25]. While shear tends to suppress the development of fingering structures that vary in the streamwise direction, it has no effect on streamwise invariant perturbations. As a result, the initial instability takes the form of sheets aligned with the flow, rather than fingers [106]. In the nonlinear regime, however, these sheets rapidly break up into 3D fingers, that are tilted in the direction of the shear. That tilt has two effects: to reduce the horizontal lengthscale of the fingers (by squeezing them together), and to add a horizontal component to the flow within the fingers. The Brown et al. [21] model for mixing by fingering convection can easily be extended to take both effects into account [25], and ultimately predicts that the mixing coefficient D fing must be corrected to account for shear in the following way: where D fing (S = 0) is the value of D fing in the absence of shear, S is the local shearing rate, λ is the growth rate of fastest-growing unsheared fingers, and χ is a number of order unity 12 . This prescription holds only for moderate shearing rates, i.e. S < 3λ or so. Beyond that, the shear itself begins to contribute to the vertical mixing as expected, but no quantitative model of the process exists (so far).    Finally, sheared fingering convection can also contribute to substantial momentum transport. We have found empirically that ν turb 0.25D fing in the limit where S < 3λ, growing to ν turb D fing for larger shearing rates. This can be used to inform models of angular momentum transport due to fingering convection in stars, which may (or may not) be useful in modeling the observed evolution of their interior rotation rates, as in [111] for example.