© ESO, 2013

1. Introduction

In models of stellar structure, situations are found where the heavier products of nuclear burning provide stability to a zone that otherwise would be unstable to convective overturning. Such a zone, or part of it, would become convective if something managed to mix its composition (Tayler 1953). The question whether such a zone should be treated as if it were mixed or not has become known as the semiconvection problem. Answers to this question differ substantially. In practice, recipes are used containing a free parameter that allows the degree of mixing to be varied. Calculations in which such a parameter is adjusted to match observations are then called “with semiconvection”. Commonly used prescriptions are those of Langer (1985) and Maeder (1997).

The presence of a semiconvective zone only has a minor effect on the thermal structure of the star. The assumed amount of mixing of composition is critical, however, because the evolution of the star is sensitive to the precise distribution of products of nuclear burning with the depth in the star. The main goal of a theory for semiconvection is thus a good determination of the rate of mixing. From the perspective of the stellar evolutionist, the theory would ideally provide a formula for the rates of mixing and transport (the effective diffusivities), as functions of local thermodynamic state and composition, and their gradients.

These formulas were derived in Spruit (1992, hereafter S92) by adapting the known physics of double-diffusive convection (Turner 1979, 1985; Proctor 1981; Huppert & Turner 1981; Schmitt 1994) to the case of a stellar interior (discussed first in this context by Spiegel 1969, 1972). The expression developed in S92 makes use of simplifications valid in the limiting case of very large Rayleigh number and very low solute diffusivity. In a companion paper (hereafter S13), the analysis in S92 is extended to cover the more moderate conditions accessible with numerical simulations. In the following, mixing rates and their dependence on the intrinsic parameters governing the problem are measured with such simulations. The predictions of S13 are tested against these results, and then used to estimate the expected mixing rates in semiconvective zones of stars.

One of the predictions in S13 is the existence of a maximum density ratio (the ratio R_ρ of stabilizing solute to destabilizing thermal stratification) for which a steady layered state is possible. In a slightly different guise, this limit also figures prominently in Proctor’s (1981) analytic analysis. In this analysis, he proved that in the limit of vanishing solute diffusivity, a layered state exists at any Rayleigh number above the critical value for convection in the absence of a stabilizing solute gradient, provided the density ratio is below this critical value. The model in S13 does not predict what happens for density ratios just above or below this maximum value. The numerical results presented in Sect. 5.2 clarify how the system behaves in this case.

The development of a linear initial gradient into the final state of overturning layers separated by diffusive steps is studied with a few examples. This transient process shows the “Kato-oscillations” expected from linear theory (see Sect. 5 for an example). It is not very relevant for astrophysical application, however, since the transition to the layered state happens on a time scale (a finite number of buoyancy oscillations) that is negligible compared with the time scales of interest, and is bound to depend on the details of the initial state. Instead, the focus here is on the transport of heat and solute in a double-diffusive staircase with layers of given thickness, as a function of the intrinsic parameters of the problem.

In Sect. 2 the known physics of double-diffusive convection of the semiconvective type is reviewed, in general and from an astrophysical perspective. Section 3 describes the transport properties of a single double-diffusive layer in terms of the model of S13. The numerical methods are given in Sect. 4. In Sect. 5 the results of our parameter study are shown and compared with the predictions of this model. Application to the case of semiconvection in a 15   M_⊙ main sequence star is discussed in Sect. 6.

2. Semiconvection and double-diffusive convection

Situations where a fluid is stabilized by the density gradient due to a dissolved heavy constituent occur in nature. An example is convection in arctic oceans (fresh meltwater cooled from above, stabilized by the salts dissolved in the sea water, see the review in Schmitt 1994). Intensively studied are East-African rift lakes (Lakes Kivu, Nyos, and Mounon, cf. Schmid et al. 2010). These are heated from below by volcanic activity, which is also a source of dissolved gases (carbon dioxide and methane, hereafter the “solute”). Their density stratification is stabilized against convection by the stable gradient resulting from the weight of the carbon dioxide. Efforts to prevent catastrophic release of carbon dioxide (Lake Nyos, e.g. Sigvaldason 1989) or to enable safe commercial exploitation of methane (Lake Kivu, Nayar 2009) have led to extensive study of the fluid flows, heat flux, and mixing rates in these natural double-diffusive systems.

The gradients in temperature and solute in the East-African lakes and the arctic are observed to be “stepped”: consisting of a stack (called “staircase”) of thin layers (decimetres to decametres). Inside a layer, overturning convection keeps the composition nearly uniform, with stable gradients in temperature and composition separating the layers. The physics involved is easily reproduced under controlled laboratory conditions (Turner & Stommel 1964; Turner 1985)¹. The layers are very long-lived: of the order of months or more in the geophysical examples mentioned, which are orders of magnitude longer than the convective turnover times inside the layers.

In the stable gradients between the overturning layers the transport of the stabilizing solute takes place by diffusion instead of convection. This strongly limits the effective transport of solute through the double-diffusive staircase. Residence times of the order of 1000 yrs are inferred for the solutes in Lake Kivu, for example (Schmid et al. 2010). This is eight orders of magnitude longer than the convective turnover times in these layers. The transport of heat is also strongly reduced and this is exploited for heat storage in solar ponds (cf. Lu & Swift 2001). Similarly low fluxes have been measured in thermohaline staircases in the arctic and antarctic oceans (e.g. Padman & Dillon 1987).

Theoretically, the observed layered nature of double-diffusive systems is now understood well. Central to this understanding is the fact that linear stability analysis does not provide relevant clues to their behaviour, because the double-diffusive case of thermal convection stabilized by a slowly-diffusing solute is subcritical. That is, a nonlinear form of the instability, in the form of a stable overturning flow, already exists below the critical temperature gradient for onset of linear instability. The linear stability condition is thus not relevant for the behaviour of the system (cf. Schladow et al. 1987).

Linear instability predicts that internal gravity modes set in above some critical value of the temperature gradient, growing in amplitude by the effect of thermal diffusion: the so-called Kato oscillations (Kato 1966). Such oscillations (cf. movie in Fig. 5) transport a negligible amount of heat or solute, compared with overturning motions of the same amplitude. For the reason alone, linear stability arguments cannot be used for estimates of the mixing rate in semiconvective zones. More important, however, is the subcritical nature of double-diffusive convection. Proctor (1981) shows analytically that, within the limit of vanishing diffusivity of the solute, the layered form of convection exists whenever the Rayleigh number exceeds the critical value for ordinary convection, irrespective of the strength of the stabilizing component. This assumption of vanishing solute diffusivity is eminently satisfied in the astrophysical case and holds reasonably well in geophysical and laboratory experiments.

In Stevenson (1979), it was assumed that the nonlinear development of the overstable oscillations would lead to saturation of the wave at a finite amplitude. This has been the rationale for some prescriptions used in astrophysics (e.g. Langer et al. 1985). The assumption of saturation at finite amplitude is appropriate for the more common supercritical forms of oscillatory instability, but not for the subcritical case, where a finite amplitude state exists for parameters where the system is still linearly stable. The assumption also conflicted with observations: the laboratory and geophysical systems all showed the same characteristic form of convection in a system of overturning layers (e.g. Turner & Stommel 1964), but none that settled into finite amplitude oscillations.

In the literature on semiconvection, it is often argued that the Prandtl number in astrophysics is much lower than in the geophysical and laboratory cases. The implied assumption that the understanding of double-diffusive convection developed in these contexts can be set aside, however, is not necessary. Proctor’s analysis, for example, is largely independent of Prandtl number. It is valid as long as Pr is not larger than of order unity, and viscosity is higher than the solute diffusivity. This is also satisfied in the astrophysical case.

In physical terms, the reason for the subcritical behaviour can be understood with an energy consideration (S92). The amount of energy it takes to overturn a layer of thickness d against a stable gradient scales as d² (as in a harmonic oscillator at the buoyancy frequency of the stratification). The expense in initial energy needed per unit mass to put the system into its finite-amplitude, layered state thus vanishes as d, down to some value where damping losses stabilize the system. A small initial perturbation or an initial Kato oscillation is sufficient to provide the energy for overturning into thin layers. Once established, this layered state is a stable form of convection. In fact, reproduction of the oscillatory phase in the laboratory requires a very careful setup of the initial stratification (Shirtcliffe 1967, see also the discussion in Huppert & Turner 1981).

This agrees with the observation in laboratory experiments and geophysical systems like Lake Kivu mentioned above, i.e. that the layering first sets in at a small thickness (cf. footnote above). In this context, the formation of layers found in the simulations by Rosenblum et al. (2011) is in line with expectations. Much less well-defined is the evolution on longer time scales, in particular the question of how and on which time scales the layer thickness evolves (see Sect. 6.6 and S13).

2.1. Double-diffusive convection in stars

Double-diffusive convection in stars has traditionally been regarded as a part of physics to be treated separately from the geophysical examples, since the numerical values of controlling parameters such as the Prandtl number are quite different. Apart from the difference in equation of state, however, the hydrodynamic equations are identical. Differences in physics that might be present between the two cases are in fact not apparent in the elementary recipes for semiconvection used in stellar evolution codes. The failure of these recipes when applied to the geophysical case is traditionally not considered an argument against their application in stars. Such recipes include (a) to assume that no mixing takes place at all (the “Ledoux” recipe); another (b) to ignore the stabilizing effect of the solute gradient (“Schwarzschild” recipe, yielding a very high mixing rate); (c) some interpolation between these recipes; (d) to assume that the amount of mixing is such that the layer becomes marginally stable to overturning; and (e) the somewhat more physically motivated oscillation-based recipe of Langer et al. above.

Define the Prandtl number Pr = ν/κ_T and the Lewis number Le = κ_S/κ_T, where κ_T is the thermal diffusivity and κ_S the diffusivity of the stabilizing “solute” (He diffusing in H, say). Because of the high thermal diffusivity mediated by photons, Pr and Le are very small numbers, some eight to ten orders of magnitude less than in geophysical cases.

Such low parameter values cannot be covered realistically in numerical simulations. If τ_c is a typical convective time scale (as estimated from the superadiabaticity and pressure scale height in the semiconvective zone), very small length scales, of the order of (τ_cκ_S)^1/2 would have to be resolved to represent the interaction between flow and diffusion. This is not computationally possible at present, even in two-dimensional simulations.

Translation from the numerically accessible parameter values to an astrophysically relevant parameter range therefore requires scaling of the results over the orders of magnitude in between. A valid extrapolation cannot be found by mere intuitive inspection of numerical results, since these are too far from the target regime. On their own, the applicability of numerical results to the astrophysical case is bound to remain diffuse, or applicability as an explicit goal has to be given up altogether.

Some physical understanding that includes this asymptotic regime is necessary. A theory that accomplishes this is given in the companion paper (S13). It is made possible by making explicit use of the observed separation into layers of convective overturning between stable diffusive interfaces.

An important difference between the astrophysical and geophysical cases is the equation of state. For the thin layering expected, sound travel times are very short compared with convective time scales. As a consequence, the fluid in a semiconvective zone behaves as if it were nearly incompressible. A Boussinesq approximation can then be used for the calculations, provided a small complication is properly taken into account. Whereas in the incompressible case the convective and diffusive heat fluxes are both governed by the temperature gradient, convection and radiative diffusion are governed by different quantities in the compressible case (entropy gradient and temperature gradient). This affects, in particular, the definition of the Nusselt number as a measure of the efficiency of heat transport. An exact translation between these cases is possible (Massaguer & Zahn 1980, see below in 6.1).

A second difference concerns the status of the heat flux in formulating the problem. The effect of semiconvective mixing on the star’s structure in evolution calculations is found to be small during the semiconvective phase itself, and is somewhat independent of the way semiconvection is approximated. (This is in part because of the limited extent of a semiconvective zone.) Its effect on the radial profile of elemental composition, however, is of lasting importance for later evolutionary stages (see Langer et al. 1989).

The consequence is that, in contrast to laboratory and geophysical situations, the heat flux F, rather than the temperature gradient in a stellar model, can be considered as known. Since the radiative contribution F_r to the heat flux is known to good approximation from the thermal structure of the star, the convective heat flux F_c = F − F_r transported by semiconvection is also known. The efficiency of convection: i.e. how close the mean thermal gradient is to the adiabatic gradient, follows from the imposed heat flux (instead of the other way around as in a laboratory experiment). We return to this distinction in Sect. 6.4.

2.2. Size of the parameter space

The fluid is described by the thermodynamic variables defining its local state (e.g. temperature and density) and the material functions (e.g. pressure, diffusivities, viscosity). In addition, the gradients of the thermodynamic variables with depth and the acceleration of gravity are relevant for the properties of the flow. Taken together, these quantities form a large parameter space, and it might be concluded that realistic numerical simulations of semiconvection would have to be done for individual zones in individual stellar models.

The equations of fluid dynamics have symmetries, however, so that the independent degrees of freedom are far fewer. They can be represented by five dimensionless parameters: a Rayleigh number Ra, the layer thickness d in units of the pressure scale height H, the Prandtl number, the Lewis number, and a density ratio R_ρ that measures the ratio of the stabilizing (solute) gradient to the destabilizing thermal gradient. The behaviour of semiconvection at some point in a star can be defined in terms of these parameters.

The Boussinesq approximation corresponds to the limit ϵ = d/H ≪ 1. The pressure scale height disappears as a parameter in this limit, and all dependence on d is subsumed in Ra. This reduces the number of parameters of the problem to four.

By a fortunate coincidence, it turns out that as long as the Prandtl number is less than unity, the results are effectively independent of Pr. This further reduces the number of independent parameters to only three. Since measurement of the mixing rate in each individual case does not require a very expensive simulation, this allows a significant volume of parameter space to be covered, and a comparison with the model predictions in S13 to be made.

3. The layered state

3.1. Layer formation, layer thickness

The formation of a layer from an initially smooth and static gradient starts with a well-known oscillatory instability, with oscillation periods close to the buoyancy period of the stratification. This is followed by nonlinear development into an overturning flow. The transition can be observed in numerical simulations (Merryfield 1995; Rosenblum et al. 2011, Sect. 5.5 below), and in very carefully designed laboratory experiments, but is not seen in geophysical cases. An unperturbed smooth initial gradient in solute and temperature is an artificial case that is unlikely to be realized in nature. In addition, the intermediate oscillatory state lasts for only a finite number of oscillation periods (5–10 in the results below), a very short time scale compared with the life time of the double-diffusive layers. In the case of Lake Kivu, for example, buoyancy periods are in the range 5–30 min, the life time for individual layers of a few months.

In laboratory and kitchen table experiments the layers are initially very thin. This is understood from the energy argument above. The layers slowly merge, either by fading of contrast between neighbouring layers, or the vertical migration of interfaces towards adjacent layers. The details of this process have not been studied much (but see McDougall 1981; Young & Rosner 2000; Ross & Lavery 2009). A plausible estimate of the rate of growth of the layer thickness can be given in terms of the effective solute diffusivity of the system, however (S13, and Sect. 6.5).

In the simulations reported here, the evolution of layer thickness is not included. The thickness (d) is therefore treated as a free parameter of the problem. It enters through the Rayleigh number (in the Boussinesq limit, see 2.2 above). In the astrophysical application of imposed heat flux, however, the resulting mixing rate is effectively independent of the layer thickness (S13 and Sect. 6.5 below). This is a fortunate circumstance, since following the evolution of a stack of layers through its merging processes would require much lengthier simulations (for an example see Young & Rosner 2000).

3.2. Structure of a layer

An individual layer in the double-diffusive staircase consists of a zone of overturning convection, separated from its neighbours by stably stratified stagnant zones. In a stagnant zone the transport both of heat and solute takes place by diffusion, so that the profiles of S and T are steep and approximately linear with depth in this zone. In the overturning zone the gradients are shallow, except in the thin boundary layers at its interface with the stagnant zone. The thickness d_s of the stagnant zone is usually only a small fraction of the layer thickness d, but can be as high as 20% close to marginal conditions for layer formation (for a discussion of its dependence on parameters see S13).

Apart from deformations by internal gravity waves in the stagnant zone, both the vertical and the horizontal velocity components vanish at its interface with the overturning zone. Apart from the small amount of solute carried by the flow, the overturning zone thus behaves essentially like laboratory convection in a box, the interfaces with the stagnant zones acting like the top and bottom plates, with viscosity enforcing no-slip conditions.Because of the periodic boundary conditions in the horizontal direction, free slip conditions would allow the stagnant zones as a whole to start moving sideways. To the overturning flow, the stagnant zones would still present internal no-slip interfaces, however. We have chosen the no-slip conditions since the freedom of such large-scale flows is somewhat of an artefact of the periodic boundary conditions in the horizontal direction. In layers of realistically large horizontal extent they would probably not happen. This picture has already been put forward early on in the interpretation of laboratory and geophysical observations (Shirtcliffe 1967; Linden & Shirtcliffe 1978) and of the results from numerical simulations, such as those reported below. It has become the standard interpretation used for deducing heat and solute fluxes in thermohaline staircases in geophysics (e.g. Padman & Dillon 1987; Turner & Stommel 1964; Schmid et al. 2010).

3.3. The overturning zone

At the boundary with the stagnant zone, the overturning flow has three nested boundary layers: for temperature, solute, and flow speed, there are the thermal, solute, and viscous boundary layers, respectively. Their thicknesses are determined by the thermal diffusivity κ_T, solute diffusivity κ_S, and (kinematic) viscosity ν. The highest diffusivity (thermal) has the widest boundary layer. For astrophysical conditions, κ_T ≫ κ_S,ν, so the solute and viscous sublayers are thinner than the thermal boundary layer.

The horizontal velocity vanishes in the stagnant zone (apart from internal wave modes), so it imposes no-slip conditions on the overturning flow. Since the viscous boundary layer is thin compared with the thermal boundary layer; however, this boundary condition has little effect on turnover times and the resulting heat flux, and it becomes noticeable only at low Rayleigh numbers.

The need to properly resolve the narrowest of the boundary layers in a numerical simulation determines the number of grid points needed in the vertical direction. In the astrophysical case, viscosity is typically similar to or somewhat larger than solute diffusivity, so the required resolution is set by the solute boundary layer (cf. Sect. 4.4).

3.4. Heat transport

The flow in the overturning zone is driven entirely by boundary layers at the top and bottom steps of the layer, much like in laboratory convection in a box (Niemela et al. 2000). Except for these thin boundary layers, the (horizontal averages of) entropy and composition are almost uniform inside the layers. Under the conditions in a stellar interior, the thermal diffusivity is much higher than the diffusivity of the solute (helium in hydrogen, say). The thickness of the solute boundary layers is then much smaller than the thermal boundary layers, and the amount of solute in transit across the layer vanishes in this limit. The convective flux is then almost the same as in the absence of a stabilizing solute (cf. Schmitt 1994). Under these asymptotic conditions, the dependence of heat flux on Rayleigh number can be taken from a simple estimate, as was done in S92, or else a result from laboratory measurements of convection can be used (e.g. Niemela et al. 2000). With the convective flow thus known, the flux of solute can be calculated as well.

Convection experiments in gaseous Helium at temperatures just above the critical point by Castaing et al. (1989) yielded the fit $\mathrm{N}{\mathrm{u}}_{\mathrm{T}}\mathrm{=}\mathrm{0.23}{\mathrm{Ra}}^{\mathrm{0.282}}$(1)over the range 10⁷ < Ra < 10¹⁴. More recent measurements by Niemela et al. (2000) using superfluid ³He, over the remarkable range 10⁶ ≤ Ra ≤ 10¹⁷ gave a marginally different result: $\mathrm{N}{\mathrm{u}}_{\mathrm{T}}\mathrm{=}\mathrm{0.124}\hspace{0.17em}{\mathrm{Ra}}^{\mathrm{0.309}}\mathit{.}$(2)This can be compared with the estimate based on a two-dimensional argument in S92: $\mathrm{N}{\mathrm{u}}_{\mathrm{T}}\mathrm{\approx }\mathrm{0.5}\hspace{0.17em}{\mathrm{Ra}}_{\mathrm{\ast }}^{\mathrm{0.25}}\mathit{.}$(3)Since the Prandtl number in the laboratory experiments (of order 0.7) is not far from unity, Ra in (1, 2) can be identified with Ra_∗ in (3). Expression (3) then agrees within a factor 2 with (2), up to Ra_∗ = 10¹⁵.

The expressions above take the form Nu_T = a   Ra^b. They are fits for large Ra; to better cover lower Rayleigh numbers as well, we modify the expression slightly: $\mathrm{N}{\mathrm{u}}_{\mathrm{T}}\mathrm{-}\mathrm{1}\mathrm{=}\mathit{a}\hspace{0.17em}\mathrm{(}\mathrm{R}{\mathrm{a}}_{\mathrm{\ast }}\mathrm{-}\mathrm{R}{\mathrm{a}}_{\mathrm{\ast }\mathrm{c}}{\mathrm{)}}^{\mathit{b}}\mathit{,}$(4)where Ra_∗ c = Pr   Ra_c, and Ra_c ≈ 1400 the critical Rayleigh number for convection with no-slip boundary conditions. With this change the Nusselt number approaches the diffusive value of 1 as Ra approaches the critical value.

3.5. Solute transport

The transport of solute is governed by the diffusion from the overturning flow into the steep gradient in the stagnant zone. This happens in a narrow boundary layer of the flow. Its thickness δ_s is of the order of δ_S = (κ_Sτ)^1/2, where τ is the time the flow travels in contact with the boundary before descending/ascending back into the interior. In the same way, the thermal boundary layer has thickness δ_T = (κ_Tτ)^1/2. The fluxes of solute and heat carried in these boundary layers also determines the flux across the overturning zone as a whole. This predicts that the solute and heat fluxes F_S, F_T are related by

${\mathit{F}}_{\mathrm{S}}\mathit{/}{\mathit{F}}_{\mathrm{T}}\mathrm{~}\mathrm{(}{\mathit{\kappa }}_{\mathrm{S}}\mathit{/}{\mathit{\kappa }}_{\mathrm{T}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}$(5)This is the well-known relation derived in various ways for double-diffusive convection (cf. Turner 1985). In a more accurate form it appears in the model of S13 and in the numerical experiments reported below.

The amount of solute transported across the overturning zone is limited by the fact that convective plumes only develop from material with a net buoyancy of the unstable sign. This limits the density contrast of the stable solute that can be carried to a fraction 1/R_ρ of the density contrast across the layer as a whole (Schmitt 1994, S92, S13). Solute with a higher density contrast is not mixed into the convective flow, but remains in the stagnant zone. The solute flux carried by convection in the overturning zone has to match the diffusive flux in the stagnant zone, and idem for the heat flux. In S13 it was shown that these conditions, together with expression (4), form a complete model, determining the effective transport coefficients, as well as the thickness of the stagnant zone, as functions of the three parameters of the Boussinesq problem, R_ρ, Ra_∗, Le.

In the estimates above the amount of solute that is transported is given by the amount flowing through the solute boundary layer of thickness δ_s. This assumes that the interaction with the stagnant zone is determined entirely by diffusion of the solute. This is a minimum: the stagnant zone as a whole is stable in the sense of the Richardson condition (one verifies that this is equivalent to the fact that R_ρ > 1), but the density gradient is lower near its boundary with the overturning zone, and a certain degree of mixing by shear instabilities is to be expected. This will increase the amount of solute that has low enough buoyancy to be carried with the flow. As argued in S13, this additional amount scales with δ_s, so that its net effect is equivalent to an increase in δ_s by a factor q, of order unity, increasing the solute flux by the same factor. The predicted relation between solute and thermal Nusselt numbers was derived in S13, and is $\mathrm{N}{\mathrm{u}}_{\mathrm{S}}\mathrm{-}\mathrm{1}\mathrm{=}\frac{\mathit{q}}{{\mathrm{Le}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{R}}_{\mathit{\rho }}}\hspace{0.17em}\mathrm{(}\mathrm{N}{\mathrm{u}}_{\mathrm{T}}\mathrm{-}\mathrm{1}\mathrm{)} \mathrm{(}{\mathit{R}}_{\mathit{\rho }}\mathit{<}{\mathrm{Le}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{)}\mathit{.}$(6)The value of q will be a fitting parameter in the quantitative comparison with the numerical results. (The value R_ρ = Le^− 1/2 is the highest for which the model predicts existence of a layered state at any Rayleigh number, see also Sect. 5.2 below.)

3.6. Model summary

The model used for comparison with the numerical simulations is defined uniquely by the ingredients described above. It makes use of only minimal assumptions about a double-diffusive layer: (a) convection in the overturning zone of a layer is described by a fit to laboratory measurements, (b) transport in the stagnant zone is by diffusion only, and (c) the fluxes of solute and heat are related by the buoyancy limit in the overturning zone (S13 Eq. (17)). These together determine the thickness of the stagnant zone and the effective transport coefficients of heat and solute through the layer as functions of the parameters of the problem, Le, Ra, and R_ρ. The behaviour of the model is discussed in more detail in S13.

4. Numerical simulations

All equations were calculated on a two-dimensional rectilinear Cartesian grid in terms of finite differences. The horizontal coordinate is labelled with x, the vertical coordinate with z, where z = 0 is the bottom and z = 1 the top boundary. The set of equations are implemented into the ANTARES software framework (Muthsam et al. 2010). Advective currents are treated by a weighted essentially nonoscillatory (WENO) finite volume scheme at fifth order (Shu & Osher 1988), and the physical diffusion is handled by a fourth-order finite difference discretization. A second-order total variation diminishing scheme is chosen as time integrator. To avoid odd-even decoupling, a MAC grid, which locates vector variables at cell faces and scalar variables at cell centres, is used.

For a detailed description of the numerical solution of the binary mixture equations presented here, see Zaussinger (2011).

4.1. Boussinesq approximation for a binary mixture

In the Boussinesq approximation the continuity equation is reduced to that of an incompressible fluid, $\frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{0}\mathit{,}$(7)such that in the equation of motion the density is taken to be a constant ρ₀, $\frac{\mathrm{d}{u}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}\mathrm{\nabla }\mathit{P}\mathrm{+}\mathrm{(}\mathrm{-}{\mathit{\alpha }}_{\mathrm{T}}\mathrm{\Theta }\mathrm{+}{\mathit{\alpha }}_{\mathrm{S}}{\mathit{S}}^{\mathrm{)}}{g}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\nu }\mathrm{\nabla }{u}\mathrm{\right)}\mathit{,}$(8)where the “expansion coefficients” α_T,α_S describe the density effects of variations Θ and S in temperature and solute, which is assumed to be low compared with the mean temperature and salinity . They are given by advection-diffusion equations: $\begin{array}{ccc}\frac{\mathrm{d\Theta }}{\mathrm{d}\mathit{t}}& \mathrm{=}& \\ \frac{\mathrm{d}\mathit{S}}{\mathrm{d}\mathit{t}}& \mathrm{=}& \end{array}$where ν is the kinematic viscosity. The mean gradients and can be expressed more usefully in terms of the buoyancy frequencies: $\begin{array}{ccc}{\mathit{N}}_{\mathrm{T}}^{\mathrm{2}}& \mathrm{=}& \\ {\mathit{N}}_{\mathrm{S}}^{\mathrm{2}}& \mathrm{=}& \end{array}$The density ratio R_ρ is then defined as in the compressible case, Eq. (32). The numerical algorithm solving this set of equations is based on a semi-implicit scheme. Intermediate values for the velocity field u^∗ are calculated explicitly from the equations of motion. By the nature of the incompressible equations, the pressure update is done implicitly, by solving a Poisson equation: $\mathrm{\Delta }\mathit{P}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{0}}}{\mathrm{\Delta }\mathit{t}}\mathrm{(}\mathrm{\nabla }\mathrm{·}{{u}}^{\mathrm{\ast }}\mathrm{)}\mathit{.}$(13)The resulting pressure P leads to the required divergence free velocity field at the new time step n + 1. ${{u}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{=}{{u}}^{\mathrm{\ast }}\mathrm{-}\frac{\mathrm{\Delta }\mathit{t}}{{\mathit{\rho }}_{\mathrm{0}}}\mathrm{\nabla }\mathit{P}\mathit{.}$(14)

4.2. Compressible fluid equations for a binary mixture

Verification of the Boussinesq results has been done with simulations of the fully explicit compressible fluid equations. The fluid is assumed to be an ideal gas, which is a good approximation to a binary gas mixture of our interest. These are the continuity equation $\frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\rho }{u}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}$(15)the partial density equation $\frac{\mathit{\partial }\mathrm{\left(}\mathit{\rho c}\mathrm{\right)}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\rho c}{u}\mathrm{\right)}\mathrm{=}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\rho }{\mathit{\kappa }}_{\mathit{c}}\mathrm{\nabla }\mathit{c}\mathrm{\right)}\mathit{,}$(16)the momentum equation $\frac{\mathit{\partial }\mathrm{\left(}\mathit{\rho }{u}\mathrm{\right)}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{(}\mathit{\rho }{uu}\mathrm{+}\mathit{P}\underline{\mathrm{I}}\mathrm{)}\mathrm{=}\mathit{\rho }{g}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathit{\tau ,}$(17)the total energy equation $\frac{\mathit{\partial \rho E}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left[}{u}\mathrm{\right(}\mathit{\rho E}\mathrm{+}\mathit{P}\mathrm{\left)}\mathrm{\right]}\mathrm{=}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}{\mathit{K}}_{\mathrm{h}}\mathrm{\nabla }\mathit{T}\mathrm{\right)}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}{u}\mathit{\tau }\mathrm{\right)}\mathrm{+}\mathit{\rho }{gu}\mathit{,}$(18)and the equation of state $\mathit{P}\mathrm{=}\frac{\mathrm{ℛ}\mathit{\rho T}}{\mathit{\mu }\mathrm{\left(}\mathit{c}\mathrm{\right)}}\mathit{,}$(19)where c is the solute mass fraction, $\mathit{E}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{|}{u}{\mathrm{|}}^{\mathrm{2}}\mathrm{+}\mathit{e}$ is the total specific energy in units of energy per mass, e the specific internal energy, K_h the heat conduction coefficient, κ_c the solute diffusion coefficient, τ the viscous stress tensor, g the gravitational acceleration, I the unit tensor, μ = μ(c) the mean molecular weight, and ℛ the gas constant.

4.3. Units, boundary and initial conditions

As unit of length we use the layer thickness d for the Boussinesq cases and the pressure scale height H for the compressible calculations. A nominal convective turnover time is used as unit of time. As units of temperature (potential temperature Θ in the compressible case) we use 1/α_T and 1/α_S for solute concentration. The density ratio then becomes . The Boussinesq equations are invariant to arbitrary additive constants in temperature and solute. We set these such that T and S are zero at the top boundary.

Because of the symmetries of the problem, there are fewer independent parameters than physical variables describing it. So some of the physical quantities appearing in the problem can be set to unity. For these we choose the temperature difference between top and bottom of the layer, the density at the bottom of the layer, and the acceleration of gravity. The Rayleigh number Ra_∗ is then controlled through the thermal diffusivity κ_T, the solute difference between top and bottom through the density ratio R_ρ, and the solute diffusivity through the Lewis number Le.

Most of the calculations were done in a box simulating a single layer from the double-diffusive staircase, so the top and bottom boundaries coincide with the steps between layers. This ignores the distortions of the interfaces by surface waves, but since the essence of the double-layering phenomenon is that the transport across the interface is by diffusion, this is not expected to make a big difference. To check that this is indeed the case, a smaller set of simulations was done in which a step is present inside the volume (Sect. 5.7). The vertical boundary conditions are thus taken to be impermeable and stress-free. In the horizontal direction periodic conditions are used: $\begin{array}{ccc}{\mathit{u}}_{\mathit{z}}& \mathrm{=}\mathrm{0} \mathrm{(}\mathit{z}\mathrm{=}& \mathrm{0}\mathit{,}\mathrm{1}\mathrm{)}\\ \frac{\mathit{\partial }{\mathit{u}}_{\mathit{x}}}{\mathit{\partial z}}& \mathrm{=}\mathrm{0} \mathrm{(}\mathit{z}\mathrm{=}& \mathrm{0}\mathit{,}\mathrm{1}\mathrm{)}\\ \mathit{S}& \mathrm{=}{\mathit{R}}_{\mathit{\rho }}\mathrm{(}\mathit{z}\mathrm{=}& \mathrm{0}\mathrm{)}\\ \mathit{T}& \mathrm{=}\mathrm{1}\mathrm{(}\mathit{z}\mathrm{=}& \mathrm{0}\mathrm{)}\\ \mathit{S}& \mathrm{=}\mathrm{0}\mathrm{(}\mathit{z}\mathrm{=}& \mathrm{1}\mathrm{)}\\ \mathit{T}& \mathrm{=}\mathrm{0}\mathrm{(}\mathit{z}\mathrm{=}& \mathrm{1}\mathrm{)}\mathit{.}\end{array}$As initial condition the stratification of temperature and solute was taken to be horizontally uniform with either a linear gradient between the values at the top and bottom (the “linear” case below), or something approximating the boundary layer structure expected of the final state (the “step” case). Small initial random perturbations are applied on the solute field.

The numerical algorithm for setting up the initial conditions for the compressible fluid equations is an extension of the procedure presented in Muthsam et al. (1995, 1999).

4.4. Numerical setup

As the thinnest boundary layer in the problem, δ_S determines the numerical resolution needed near the boundaries. Since the fine structure in the interior of the layer largely consists of boundary layers “peeled off” from the boundaries, the same resolution is needed in the interior of the layer as well, and there is no need or justification for using non-uniform grids.

The expected solute boundary thickness for Ra_∗ = 1.6 × 10⁵ and a Lewis number of Le = 0.1 is 1.6% of the layer thickness d. With the high-order spatial discretization used, a nominal resolution of n_BL = 3 points is needed across the critical structures that must be captured. The transport of solute across the layer is determined by the vertical structure of the boundary layers. Resolving these with n_BL = 3 translates to the n_z = 200 points needed in the vertical direction for this case.

To test the actual convergence of the results with respect to resolution, a series of simulations with n_z = 100, 200, 300 and 600 has been run for Ra_∗ = 1.6 × 10⁵, Le = 0.1, and R_ρ = 1.15, see Table 1. Convergence to 15% is already reached at n_z = 100. For the lowest Lewis number used in the results reported below in Table 2, Le = 0.01, the number of grid points needed would be $\sqrt{\mathrm{1}}\mathrm{0}$ times higher. The number of points across the height of the box used in these simulations, n_z = 300, is deemed sufficient for the level of accuracy we are aiming for.

Most of the calculations were done with a horizontal-to-vertical aspect ratio of 2:1. A few tests with different ratios showed that this choice does not affect the measured Nusselt numbers significantly (Sect. 5.4).

5. Results of numerical simulations

Over a wide range in parameter values, about 100 numerical simulations were done in the Boussinesq approximation and about 20 simulations were performed with the fully compressible code. Only the cases closest to the parameter range of astrophysical interest (about 60) are reported here (Figs. 2, 4, Table 1).

An example of the flow structure is shown in Fig. 1, with Pr = 0.1, Ra_∗ = 5 × 10⁵, Le = 0.01, and R_ρ = 1.15. It shows the key characteristics of double-diffusive convection: the boundary layers and the plumes of the solute are narrower than those of temperature, on account of the low solute diffusivity.

Fig. 1 Flow structure in a double-diffusive layer. The temperature field (top) is more diffuse than the solute (“helium”) concentration, as a result of the high thermal diffusivity.

5.1. Dependence on RA_∗ and R_ρ

Figure 2 shows the dependence of the measured Nusselt numbers Nu_T and Nu_S on the parameters Ra_∗ and R_ρ, for the case Le = 0.01 and Pr = 0.1. The results show some fluctuation, as expected from the limited number of overturning times for which the simulations were run. This results in uncertainties in the Nusselt numbers of the a factor 1.5. For comparison with the theoretical model predictions (dashed lines) the adjustable erosion factor q in Eq. (6) has been set to 1.5. Within the scatter, the predictions appear consistent with the numerical results to a factor of 2 or less. In fact, leaving out the erosion factor (i.e. setting q = 1) does not make the fit much worse.

Fig. 2 Dependence of the Nusselt numbers on density ratio Rρ and Ra∗, for Le = 0.01 and Pr = 0.1. Solid: numerical results; dashed: model predictions for erosion factor q = 1.5.

The upper right-hand panel shows Nu_s vs. Nu_T, where each curve represents the results for a fixed density ratio R_ρ and concurrently increasing Ra_∗. The approximately linear relation between the two shows that the ratio of thermal-to-solute transport depends primarily on Lewis number and density ratio, but not much on Ra, as predicted by (6). The actual transport efficiency of both depends, of course, on the Rayleigh number Ra_∗ as well (upper left panel).

5.2. Behaviour near the maximum density ratio

The data in Fig. 2 show how the Nusselt number Nu_T declines with increasing density ratio. The theoretical analysis in S13 predicts a maximum density ratio R_ρ   max for the existence of a (statistically) stationary layered state. It is a slowly increasing function of Ra_∗, asymptotically approaching the value Le^− 1/2. To test this a series of simulations was done at Ra_∗ = 10⁶, Le = 0.1. The predicted value of R_ρ   max for this combination is 1.2 (the asymptotic value would be 1.3). Figure 3 shows the resulting Nusselt numbers as functions of time.

For R_ρ = 1.2 and 1.3, the heat flux in the simulations quickly settles to a statistically steady value after a transient due to the initial conditions used. For the three highest values of R_ρ, the initial state smoothly settles in to the purely diffusive state N_T = 1; the faster and smoother it is, the higher the density ratio. Values in between 1.5 and 1.8 are characterized by large fluctuations superposed on a slow general decline, as if the system vacillates between an overturning or a static. A precise boundary for deciding when an overturning state can exist, is therefore somewhat hard to define: it appears to depend on the length of time over which the flow is followed.

The time dependences of the Nusselt number shown in Fig. 3 suggest a value R_ρ   max ≈ 1.4, somewhat higher than the theoretical value of 1.2. For density ratios slightly higher than 1.2, it looks at first as if a statistically steady state has been reached, but on a longer time scale, a continuing decline of Nu_T is observed. The behaviour with increasing R_ρ is therefore somewhat gradual, at least for the finite length of time over which the simulations have followed the development of the flow. Interestingly, the distinction between diffusive and overturning final states becomes sharper with increasing length of time over which the simulation is followed. This also affects the results in Fig. 2. At the highest density ratios shown there, the numerically measured Nusselt numbers are systematically higher than predicted. This can be attributed to the fact that at these R_ρ the simulations were not run as long as those of Fig. 3.

The numerical results thus confirm the predicted maximum density ratio for the overturning layered state, but with the added twist of very long settling times when R_ρ is near the maximum value. This also answers a question that was left undecided in S13: the theory does not predict what happens close to the maximum density ratio. The results in Fig. 3 indicate that the system vacillates between overturning and diffusive states for increasingly long periods of time as the maximum is approached.

The fast settlement towards a diffusive state for high values of R_ρ is not surprising given that the linear stability criterion predicts exactly such a result for R_ρ > (1 + Pr)/(Le + Pr) (e.g Veronis 1968; Stevenson 1979). For the parameter values of Fig. 3, this yields R_ρ ≳ 1.8. The difference between 1.8, where Nu_T − 1 drops only initially, but then only mildly so if at all, and a value of 2.0, which shows a continuous decay seems consistent with this classical result.

The results also provide an interesting link to Proctor’s (1981) analysis. In this theory the maximum value R_ρ = Le^− 1/2 appears as a limit for its applicability. Proctor’s analysis therefore does not answer what happens near this critical number, any more than the model in S13 does, but it is gratifying that it appears to have a consistent place in both, as well as in the numerical results.

Fig. 3 Time dependence of the thermal Nusselt number NuT in simulations with Pr = 1, Le = 0.1, and Ra∗ = 106, for density ratios (labeled) near the maximum value for which a layered state is predicted by the theoretical model. Time t in units of the thermal buoyancy time scale ${\mathit{N}}_{\mathrm{T}}^{-1}$.

5.3. Dependence on Pr

Figure 4 shows the results of a set of simulations testing the dependence on Prandtl number. The Lewis number is fixed at Le = 0.01. The density ratio is varied over the range 1.15 < R_ρ < 5, Prandtl number from Pr = 0.01 to Pr = 1. Within the limited numerical sensitivity due to the stochastic nature of the flow, little systematic dependence on Pr is detectable. As argued above, a low dependence on Pr was to be expected in the limit Pr ↓ 0. It appears that this also extends to a Prandtl number of order unity.

Fig. 4 Dependence on Prandtl number and density ratio. Solute Nusselt number (top) and thermal Nusselt number (bottom) for Pr = 1, 10-1, and 10-2, with Ra∗ = 1.6 × 105 and Le fixed at 0.01.

5.4. Dependence on aspect ratio

The aspect ratio of 2:1 used in the results reported above was tested against 5:1 and 10:1 at the same spatial resolution. For the reference simulation (Pr = 0.1, Le = 0.01, Ra_∗ = 10⁵, R_ρ = 1.15), we find Nu_S = 90 and Nu_T = 8.5. By comparison the simulation with 5:1 results in Nu_S = 90 and Nu_T = 9.0. The most extended box with an aspect ratio of 10:1 and a spatial resolution of 1500 × 300 has Nusselt numbers of Nu_S = 80 and Nu_T = 8.75. The aspect ratio thus has no significant influence on the dependences of the fluxes on input parameters in our simulations, within the fluctuations due to the stochastic nature of the flow.

5.5. Dependence on initial stratification

As initial state we used either a constant linear gradient of temperature and solute between top and bottom values (“linear”), or a profile approximating the expected steplike combination of a stagnant and an overturning zone. The linear case shows how the oscillatory phase due to the Kato instability develops into overturning flow, see Figs. 5 and 6. The duration of the initial formation process is mainly determined by R_ρ and Ra_∗. It is in the range from five to ten oscillation periods of the initial stratification. The end state in both cases is statistically the same. The “step” as initial condition also covers cases that are stable in linear theory because of the subcriticality of the system, and would not develop from a linear initial profile. With this initial state, computing time can be saved for low values in Ra_∗ and high R_ρ. It was used in most of the results reported above.

Fig. 5 Example of the development of an overturning flow from Kato oscillations for Pr = 1, Le = 10-1, Rρ = 1.15, and Ra∗ = 1.6 × 105 starting from a linear stratification. Time from left to right and top to bottom. Wave braking occurs after 5 oscillations (Fig. 5.2). The layer is fully evolved after 10 oscillations (Fig. 5.4). See also the movie online.

Fig. 6 Top panel: evolution of the Nusselt number for the linear initial stratification. Convective cells get established from Kato oscillations (cf. Fig. 5) after about 10 turnover times. Starting the simulation from a step (bottom) saves computing time. At the end of the runs the Nusselt numbers of both simulations are the same within the statistical variations.

5.6. Comparison with compressible results

The compressible simulations are based on a fifth-order weighted ENO scheme. Compressible fluids lead to restrictions in time stepping (due to the need to resolve sound waves). The compressible simulations take up to 100 times longer for the same resolution compared to the incompressible solver. In both regimes a resolution of 300 × 300 points has been set. Therefore only a few tests have been done with the compressible code. The degree of compressibility is governed by the ϵ = d/H of layer thickness to pressure scale height. In the limit ϵ → 0, there is a direct translation between the compressible and the Boussinesq case (Sect. 6.1). The results of a numerical comparison with ϵ = 0.1 is shown in Table 1. The resulting Nusselt numbers do not differ significantly.

Simulations done with ϵ = 1.0 behave quite similar to those done with ϵ = 0.1. The mixing processes do not significantly differ, at least for Pr ≤ 0.1 and as long as the Rayleigh number is high enough, Ra > 5.0 × 10⁵. Differences of around 50% in the Nusselt numbers were found, but considering the present level of accuracy (a factor of 2, say), we conclude that the Boussinesq approximation qualitatively gives the right results, even for layer thicknesses approaching a scale height. At Pr of order unity or larger, the effect of viscous damping might become more significant in more strongly stratified cases, however.

5.7. Multilayer simulations

In all of the above we have assumed that the interfaces between the double-diffusive layers can be approximated as solid boundaries. To test the reliability of this assumption, a few cases were run where the initial state consisted of two steps instead of a single one. In some, though not all of these runs, the division into two layers remained until the end. An example is shown in Fig. 7 for a case with Pr = 1, Le = 10^-2, R_ρ = 1.15, and Ra_∗ = 6 × 10⁵ and a resolution of 500 × 500. We note the approximate (anti-)symmetry of the plumes near the interface in the middle, a phenomenon known from laboratory experiments (e.g. Fig. 1 in Turner 1985). It is caused by the continuity of the horizontal velocity across the interface enforced by viscosity. This symmetry is less marked in simulations at lower Prandtl numbers, when the thickness of the viscous boundary layer becomes smaller than the thickness of the stagnant zone.

Figure 8 shows horizontal and temporal average profiles of temperature and solute with height. The transition in the middle is broad compared with the boundary layers at top and bottom. Inspection of the time-dependent flow shows that this is due to two separate effects. One is the displacements of the interface by surface waves, which smoothen the average gradient without changing the actual fluxes across the interface. In addition there is possibly some real mixing associated with breaking of the surface waves, but it remains localized around the interface between overturning and stagnant zone.

Fig. 7 Snapshot of a simulation including the free interface between two layers. Left: solute, Right: temperature. Pr = 1.0, Le = 0.01, Rρ = 1.15, Ra∗ = 6 × 105. See also the movie online.

Fig. 8 Average temperature and solute profiles with height z in the double-layer simulation of Fig. 7.

6. Application to stars

6.1. Boussinesq limit: thin layers

If the double-diffusive layering is thin, convective overturning times are longer than the sound crossing time of the layer thickness, and a Boussinesq (or other “low Mach number” approximation) can be used for a compressible fluid. By the presence of a stratification in density and temperature, however, a heat flux is present even in a convectively neutral stratification. This requires some care in translating Boussinesq results to the astrophysical case.

6.2. Heat flux

We let F_r,c be the radiative and hydrodynamic (semiconvective) contributions to the heat flux in the star. The radiative heat flux is proportional to the logarithmic temperature gradient ∇ = dlnT/dlnP: ${\mathit{F}}_{\mathrm{r}}\mathrm{=}\mathit{k}\mathrm{\nabla }\mathrm{=}\mathit{k}{\mathrm{\nabla }}_{\mathrm{a}}\mathrm{+}\mathit{k}\mathrm{(}\mathrm{\nabla }\mathrm{-}{\mathrm{\nabla }}_{\mathrm{a}}\mathrm{)}\mathrm{\equiv }{\mathit{F}}_{\mathrm{ra}}\mathrm{+}{\mathit{F}}_{\mathrm{rs}}\mathit{,}$(26)where k is a constant depending on the local thermodynamic state, ∇_a the adiabatic gradient, and F_ra, F_rs are the contributions to the radiative heat flux of the adiabatic and the superadiabatic parts of the mean temperature gradient, respectively. In the Boussinesq model, the contribution F_ra is absent: the convective and radiative heat fluxes F_c and F_r are governed by the same temperature gradient. Related to this, the Boussinesq model has one parameter less: the pressure scale height H. Because of this difference, the heat flux in the stellar (compressible) model cannot be compared directly with the heat flux in an incompressible model. Instead, the ratio of convective to radiative heat flux in the Boussinesq model is to be identified with the ratio f = F_c/F_rs of the convective flux to the superadiabatic component of the radiative flux in the star, in the limit H → ∞ (cf. Massaguer & Zahn 1980). The Boussinesq model is thus the limit ϵ ↓ 0 (taken at fixed Ra_∗). If the semiconvective layer thickness ϵ is sufficiently small, the compressible case can be compared directly with the Boussinesq model (as verified by the numerical tests with ϵ = 1 and ϵ = 0.1 in Sect. 5.6 above).

This makes the semiconvection problem more amenable to numerical simulation. In contrast to a convective stellar envelope, for example, with its many scale heights to be covered, the layered nature of double-diffusive convection puts it in a parameter range that is much closer to conditions accessible with realistic ab initio calculations.

With (26) for the radiative flux, the total heat flux can be written as $\mathit{F}\mathrm{=}{\mathit{F}}_{\mathrm{r}}\mathrm{+}{\mathit{F}}_{\mathrm{conv}}\mathrm{=}\mathit{k}{\mathrm{\nabla }}_{\mathrm{a}}\mathrm{+}\mathrm{Nu}\hspace{0.17em}\mathit{k}\mathrm{(}\mathrm{\nabla }\mathrm{-}{\mathrm{\nabla }}_{\mathrm{a}}\mathrm{)}\mathit{,}$(27)where Nu is the Nusselt number. By taking out the radiative flux associated with the adiabatic gradient, the Nusselt number as defined in (27) can now be identified with the Boussinesq equivalent (see also Massaguer & Zahn 1980). The relevant Nusselt number is thus not simply the ratio of total to radiative heat flux².

6.3. Rayleigh number and density ratio

The modified Rayleigh number for a layer of thickness d is, in astrophysical notation: $\mathrm{R}{\mathrm{a}}_{\mathrm{\ast }}\mathrm{=}\mathrm{Pr}\hspace{0.17em}\mathrm{Ra}\mathrm{=}\frac{\mathit{g}{\mathit{d}}^{\mathrm{4}}}{{\mathit{\kappa }}^{\mathrm{2}}\mathit{H}}\mathrm{(}\mathrm{\nabla }\mathrm{-}{\mathrm{\nabla }}_{\mathrm{a}}\mathrm{)}\mathit{,}$(28)where κ is the thermal diffusivity (cm²/s) and g the acceleration of gravity. Ra_∗ is typically a very large number unless the layer thickness d is small (for estimates see the 15 M_⊙ example below). Apart from Le, the problem depends on the relative strength of the stabilizing solute gradient relative to the destabilizing thermal gradient. This can be measured in terms of the thermal and solute buoyancy frequencies N_T, N_S: $\begin{array}{ccc}{\mathit{N}}_{\mathrm{T}}^{\mathrm{2}}& \mathrm{=}& \frac{\mathit{g}}{\mathit{H}}\mathrm{(}{\mathrm{\nabla }}_{\mathrm{a}}\mathrm{-}\mathrm{\nabla }\mathrm{)}\mathit{,}\\ {\mathit{N}}_{\mathrm{S}}^{\mathrm{2}}& \mathrm{=}& \frac{\mathit{g}}{\mathit{H}}\hspace{0.17em}{\mathrm{\nabla }}_{\mathit{\mu }}\mathit{,}\end{array}$with ${\mathrm{\nabla }}_{\mathit{\mu }}\mathrm{\equiv }\mathrm{d}\ln \mathit{\mu }\mathit{/}\mathrm{d}\ln \mathit{P,}$(31)where μ is the mean atomic weight per particle, and the gradients are understood as mean values over a double-diffusive step. The density ratio R_ρ, ${\mathit{R}}_{\mathit{\rho }}\mathrm{\equiv }\mathrm{-}{\mathit{N}}_{\mathrm{S}}^{\mathrm{2}}\mathit{/}{\mathit{N}}_{\mathrm{T}}^{\mathrm{2}}\mathrm{=}\frac{{\mathrm{\nabla }}_{\mathit{\mu }}}{\mathrm{\nabla }\mathrm{-}{\mathrm{\nabla }}_{\mathrm{a}}}\mathit{,}$(32)is then the equivalent of the density ratio in the Boussinesq formulation, and is the dimensionless measure we use for the relative strength of the stabilizing helium gradient. Semiconvection, i.e. a stratification between Schwarzschild and Ledoux then corresponds to R_ρ > 1. ∇ − ∇_a is typically a small number in a stellar interior (even in the presence of the double-diffusive steps), and R_ρ a large number (see below).

6.4. Extrapolation

The connection between the numerical results and the astrophysical conditions requires extrapolation that cannot at present be covered by the simulations. It is covered by the model in S13, as validated by the comparison with our numerical results. In contrast with the laboratory setup where the temperature difference is the imposed quantity, the heat flux is the fixed quantity under the conditions in a stellar semiconvective zone, since the structure of the star only depends a little on the behaviour of the semiconvective zone. The asymptotic conditions for large Ra_∗ then lead (S13 Sect. 6) to simple expressions for the superadiabaticity ∇ − ∇_a and the effective Helium diffusivity κ_s   eff: $\begin{array}{ccc}\mathrm{\nabla }\mathrm{-}{\mathrm{\nabla }}_{\mathrm{a}}& \mathrm{\approx }& \\ {\mathit{\kappa }}_{\mathrm{s}\hspace{0.17em}\mathrm{eff}}& \mathrm{=}& \end{array}$(S13 Eqs. (57), (59)). The analysis in S13 does not take a modification due to the effect of radiation pressure discussed in S92 into account. With this modification, Eq. (34) becomes ${\mathit{\kappa }}_{\mathrm{s}\hspace{0.17em}\mathrm{eff}}\mathrm{=}\mathrm{(}{\mathit{\kappa }}_{\mathrm{s}}\mathit{\kappa }{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\left(\frac{\mathrm{4}}{\mathit{\beta }}\mathrm{-}\mathrm{3}\right)\frac{\mathrm{(}{\mathrm{\nabla }}_{\mathrm{r}}\mathrm{-}{\mathrm{\nabla }}_{\mathrm{a}}\mathrm{)}}{{\mathrm{\nabla }}_{\mathit{\mu }}}\mathit{,}$(35)where β = P_g/(P_g + P_r) is the ratio of gas to total pressure. The density ratio approaches the value ${\mathit{R}}_{\mathit{\rho }}\mathrm{\approx }{\mathrm{Le}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(36)while the relative thickness of the stagnant zones in this limit is of order (setting β = 1) $\mathit{\delta }\mathrm{\approx }\mathrm{1}\mathit{/}\mathrm{N}{\mathrm{u}}_{\mathrm{s}}\mathrm{=}{\mathrm{Le}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathrm{\nabla }}_{\mathit{\mu }}\mathit{/}\mathrm{(}{\mathrm{\nabla }}_{\mathrm{r}}\mathrm{-}{\mathrm{\nabla }}_{\mathrm{a}}\mathrm{)}\mathit{.}$(37)The range of validity of the limit is ${\mathit{l}}_{\mathrm{0}}\mathrm{\ll }\mathit{d}\mathrm{\ll }\mathit{H,}$(38)where ${\mathit{l}}_{\mathrm{0}}\mathrm{=}\mathrm{(}{\mathit{\kappa }}^{\mathrm{2}}\mathit{H}\mathit{/}\mathit{g}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{4}}$(39)is the length scale on which the thermal diffusion time scale equals the free fall time over a pressure scale height.

The asymptotic value R_ρ = Le^− 1/2 (Eq. (36)) is also the maximum density ratio for which the theory predicts the existence of the layered state. Semiconvection in a stellar model with heat flux fixed is thus predicted to settle close to this maximum, if the layer thickness is not small (d ≫ d₀ or Ra_∗ ≫ 1). How close it actually settles can also be predicted, but requires analysis to the next order in the small quantity 1/Ra_∗.

In a stellar evolution model, expressions (33, 35) give nearly identical results to the model of Zaussinger (2011), which is based on S92. They are easier to implement since they have been restricted to conditions d ≫ l₀ that are the most interesting for stellar evolution anyway.

6.5. A 15 M_⊙ star

We are now in the position to estimate the range of parameter values for semiconvection in a star. Consider the important case of massive stars around main sequence turnoff. We use the model kindly provided by A. Weiss (model “fzm15_151”). Characteristic values for the physical quantities in the semiconvective zone of this model are g ≈ 10⁶ cm/s², κ_S ≈ 1 cm²/s, κ ≈ 3 × 10⁸ cm²/s, H ≈ 2 × 10¹⁰ cm, ∇_a = 0.4, ∇_r − ∇_a ≈ 0.02, and ∇_μ ≈ 1.

With Eqs. (28) and (33), the Rayleigh number can be found for this semiconvective zone, as a function of the layer thickness d. This yields Ra_∗ ~ 10¹² for d/H = 0.1, or Ra_∗ ~ 10⁸ for d/H = 0.01, for example. For such layer thicknesses, conditions are thus in the asymptotic regime for which the expressions in 6.4 are approximately valid.

The effective He-diffusivity from (35) is κ_s   eff ≈ 10³ cm²/s, three orders of magnitude above the microscopic value. The mixing time scale over a pressure scale height is thus about 10¹⁰ yr. The value of l₀ is 2    ×    10⁵ cm, or about 10^-5 pressure scale heights. The limiting expressions (33−37) are thus applicable over a large range in the value of the (uncertain) layer thickness. The lower end of this range is not likely to be very relevant, since at d ~ l₀ the expected time scale for layer merging (see below) is very short, of the order of Le^−1/2 times the free fall time scale over a pressure scale height (from Eqs. (39), (42)).

The time scale for the initial layer formation from a Kato oscillation is similarly short, a few times the instability growth time, which is close to the convective overturning time in the absence of the stabilizing μ-gradient. For the 15 M_⊙ star, this works out to about one day. The asymptotic expressions (33) and (35) should therefore be adequate for use in stellar evolution codes.

6.6. Layer thickness

While the layer thickness has little influence on the effective mixing rate, it is useful to check whether it could become significant compared with the pressure scale height, in which case the limit of thin layering assumed would not be valid. This requires a model for the layer thickness, for which no good theory is available. Observations in laboratory experiments and in geophysical cases show that layer thickness is not constant, but grows in time by processes in the merging of neighbouring layers (McDougall 1981; Ross & Lavery 2009, see also the numerical experiment in Young & Rosner 2000). Layer thickness can therefore not be treated independently of the history of the system.

An estimate can be made, however, from the effective solute diffusivity. Changes in layer thickness involve the exchange of solute between neighbouring layers. Over a distance D, the time scale τ for exchange is given by the effective diffusivity κ_s   eff: τ ≈ D²/κ_s   eff. Setting D equal to the layer thickness d itself then gives an estimate of the rate of change of thickness by merging $\mathrm{d}\ln \mathit{d}\mathit{/}\mathrm{d}\mathit{t}\mathrm{=}\mathrm{1}\mathit{/}\mathit{\tau ,}$(40)or $\frac{\mathrm{d}}{\mathrm{d}\mathit{t}}\mathit{d}\mathrm{=}\frac{{\mathit{\kappa }}_{\mathrm{s}\hspace{0.17em}\mathrm{eff}}}{\mathit{d}}\mathrm{·}$(41)If κ_s   eff is constant in time, then $\mathit{d}\mathrm{=}\mathrm{(}\mathrm{2}{\mathit{\kappa }}_{\mathrm{s}\hspace{0.17em}\mathrm{eff}}\mathit{t}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(42)where t is time since the semiconvective condition started. In the 15   M_⊙ star the duration of the semiconvective state is of the order of 10⁷ yr, so the layer thickness to be expected at the end of the semiconvective phase is of the order of 10⁹ cm, or 0.05   H. Given the uncertainties involved, it cannot be ruled out that semiconvective layer thickness d can actually become macroscopic, of order H, in the course of the star’s evolution.

7. Discussion

As in the case of ordinary convection, there are far fewer intrinsic parameters in semiconvection than the quantities defining the physical state in a stellar interior. This allows a significant volume of astrophysically useful parameter space to be covered by a grid of numerical simulations. A physical model as applied here is still needed, however, to interpret the results and to extrapolate them meaningfully to the astrophysically relevant range.

The low value of the viscosity and solute diffusivity compared with the thermal diffusivity constitute a limiting case that actually simplifies the double-diffusive problem greatly. Among other things, the results become nearly independent of the Prandtl number in this limit (as suggested already by Proctor’s 1981 analysis).

The simulations were all done in 2D, so tests in 3D will be needed for verification. It is unlikely, however, that the results will turn out very different, at least within an astrophysical factor of 2. The reason is that in the low-Pr, low-Le limit the flow in the layers is almost equivalent to ordinary convection between plates. Owing the low solute diffusivity, the amount of solute in the bulk of the layer is small. It can then be treated as a perturbation, as assumed in S92. Known laboratory results for convection at very high Rayleigh numbers can then be used to extrapolate the numerical results. As shown in S13, this makes predictions similar to the simple 2D model used in S92. In particular the effective mixing rate is very low, as also found (in a different region in parameter space) in geophysical examples, such as Lake Kivu and the (Ant)Arctic oceans (cf. Sect. 2 and S13).

Stochastic fluctuations in the flow produce scatter in the fluxes of heat and solute, which affect the measured averages (cf. Table 1). The accuracy of these averages could be improved with longer runs.

Most of the results presented are based on simulation of a single double-diffusive layer. As we found in Sect. 5.7, this does not fully reproduce broadening of the interfaces between layers by surface waves (see movie in Fig. 7).

The physics determining the thickness of the individual double-diffusive layers remains uncertain. In the equivalent geophysical examples semiconvective zones always consist of a large number of very long-lived layers. Observations in the East-African rift lakes (e.g. Schmid et al. 2010) show that layers first forming at the boundary of an expanding double-diffusive zone are always thin, subsequently growing slowly by a process of merging. This is understood in terms of an energy argument (cf. Sect. 2), and is likely to happen in a growing stellar semiconvective zone as well.

The estimate (Eq. (42)) suggests that layer thickness might approach a significant fraction of a pressure scale height. In this case, only a few layers would be left at the end of the semiconvective phase, and the location of their boundaries may well vary somewhat randomly between stars. As surmised in Zaussinger (2011), this might introduce a random element in the late stages of the evolution of massive stars.

Semiconvection as studied here addresses only one of two astrophysical meanings of the term. In addition to the effect of a stabilizing solute gradient, there is the effect of a composition-dependent opacity: the increase in the radiative gradient when mixing is assumed, with consequent onset of convection in a stratification that was radiative before mixing. Historically, this has actually been the main concern in stellar evolution computations. The physics of this kind of semiconvection is completely different from the double-diffusive kind (cf. Kippenhahn et al. 2013); it has not received the same theoretical attention.

8. Conclusions

The results show how the physics of double-diffusive convection, which is known from geophysical and laboratory studies, can be applied to the astrophysical case of low Prandtl and Lewis numbers. As expected, the process takes place in the form of the characteristic double-diffusive layering known and theoretically understood since the 1980s.

In the asymptotic regime occupied by astrophysical conditions, the number of independent parameters determining the physics is effectively reduced to three, so that a meaningful range of parameter values can be covered with numerical simulations. We have compared a grid of simulations with the predictions of the model for a semiconvective layer in S13. In the parameter region where the estimate overlaps with the region covered by the simulations we found good agreement, even without significant tuning of the one fitting parameter used (the “erosion factor” q, Sect. 3.5).

The simulations also provide an answer to an interesting question raised by Proctor’s (1981) mathematical analysis and by the model in S13. Both predicted that a critical number R_ρ = Le^− 1/2 plays a role in the steady layered state. The results presented in Sect. 5.2 clarify the behaviour of the layered state around this maximum density ratio. In a system set up in a layered state just above R_ρ   max, the Nusselt numbers show large fluctuations in amplitude, and the duration of the fluctuations lengthens as R_ρ approaches R_ρ   max. This explains why theories based on the assumption of stationarity fail near this limit.

Boussinesq and fully compressible simulations give equivalent results in the limit of small layer thickness. The Boussinesq calculations even reproduce the approximate Nusselt numbers for layers as thick as a scale height. The numerical results confirm the theoretical expectation that the dependence on viscosity is weak at low Prandtl number and extends to Pr of order unity.

The effective mixing rate in the semiconvective zone of a 15   M_⊙ main sequence star is predicted to be low, though the semiconvective phase may possibly leave significant jumps in the profile of helium concentration.

Online material

Movie to Figure 5 Access here

Movie to Figure 7 Access here

Acknowledgments

We thank the referee for a perceptive and detailed review, which has led to substantial improvements. F. Zaussinger is grateful to F. Kupka and H.J. Muthsam for discussions of the numerical solution of binary mixture equations with high resolution methods. He was supported by the DFG within the project “Modelling of diffusive and double-diffusive convection” (projects KU 1954-3/1 and KU 1954-3/2 in SPP 1276/1 and SPP 1276/2, project leader F. Kupka) within the interdisciplinary Metstroem project. He was also supported by the Austrian Science Foundation (project P20973, Numerical Modelling of Semiconvection, project leader H.J. Muthsam). We also thank H. Grimm-Strele for implementing a parallel Poisson solver in the ANTARES code suite.

References

All Tables

All Figures

	Fig. 1 Flow structure in a double-diffusive layer. The temperature field (top) is more diffuse than the solute (“helium”) concentration, as a result of the high thermal diffusivity.
In the text

	Fig. 2 Dependence of the Nusselt numbers on density ratio Rρ and Ra∗, for Le = 0.01 and Pr = 0.1. Solid: numerical results; dashed: model predictions for erosion factor q = 1.5.
In the text

	Fig. 3 Time dependence of the thermal Nusselt number NuT in simulations with Pr = 1, Le = 0.1, and Ra∗ = 106, for density ratios (labeled) near the maximum value for which a layered state is predicted by the theoretical model. Time t in units of the thermal buoyancy time scale ${\mathit{N}}_{\mathrm{T}}^{-1}$.
In the text

	Fig. 4 Dependence on Prandtl number and density ratio. Solute Nusselt number (top) and thermal Nusselt number (bottom) for Pr = 1, 10-1, and 10-2, with Ra∗ = 1.6 × 105 and Le fixed at 0.01.
In the text

	Fig. 5 Example of the development of an overturning flow from Kato oscillations for Pr = 1, Le = 10-1, Rρ = 1.15, and Ra∗ = 1.6 × 105 starting from a linear stratification. Time from left to right and top to bottom. Wave braking occurs after 5 oscillations (Fig. 5.2). The layer is fully evolved after 10 oscillations (Fig. 5.4). See also the movie online.
In the text

	Fig. 6 Top panel: evolution of the Nusselt number for the linear initial stratification. Convective cells get established from Kato oscillations (cf. Fig. 5) after about 10 turnover times. Starting the simulation from a step (bottom) saves computing time. At the end of the runs the Nusselt numbers of both simulations are the same within the statistical variations.
In the text

	Fig. 7 Snapshot of a simulation including the free interface between two layers. Left: solute, Right: temperature. Pr = 1.0, Le = 0.01, Rρ = 1.15, Ra∗ = 6 × 105. See also the movie online.
In the text

	Fig. 8 Average temperature and solute profiles with height z in the double-layer simulation of Fig. 7.
In the text