CHEMICAL TRANSPORT AND SPONTANEOUS LAYER FORMATION IN FINGERING CONVECTION IN ASTROPHYSICS

In Figure 1, we present the temporal evolution of the thermal and compositional Nusselt numbers in a typical numerical simulation of fingering convection, taking Pr = 1/10, τ = 1/30, and R₀ = 3. This evolution begins with a period of exponential growth from t = 0 to t = 160. The observed growth rate is found to be twice the growth rate of the fastest growing mode according to a linear stability analysis of the governing equations. This is as expected, as shown in Section 4.1. The transport peaks around t = 165, after which the fluxes of temperature and composition saturate (at a level slightly below the peak) and remain roughly constant for the remainder of the simulation.

Zoom In Zoom Out Reset image size

In order to visualize the saturation process, we present snapshots of the compositional field of the simulation at three stages before and after saturation in Figures 2(a)–(c), marked with vertical lines in Figure 1. In Figure 2(a), which shows the state just as the primary fingering instability begins to develop, we see tall and thin vertical structures, called "elevator modes," with high μ plumes sinking and low μ plumes rising. In Figure 2(b), which is taken just before the peak in Figure 1, we recognize the same vertical structures but note that smaller-scale instabilities have distorted them. The latter eventually destroy the elevator modes and the system achieves saturation in Figure 2(c). At this point, very little of the original elevator modes remain discernible, and quasi-steady, homogeneous turbulence dominates the dynamics.

Zoom In Zoom Out Reset image size

While most simulations behave exactly as Figure 1, in some cases where R₀ is chosen to be very close to the onset of overturning convection (i.e., for small r, see Equation (11)), the Nusselt numbers begin to increase gradually again after saturation. This behavior appears clearly in Figure 3, for which Pr = 1/10, τ = 1/30, and R₀ = 1.1. We distinguish two phases post-saturation: from t = 100 to t = 600, the mean Nusselt numbers increase slowly and undergo quasi-periodic oscillations. We find that these oscillations have twice the buoyancy frequency, which is what we expect for gravity waves.⁶ Since linear gravity waves do not lead to an increase in the mean flux, the gradual rise of the mean Nusselt numbers indicates that the waves are highly nonlinear. At t = 600, the transport suddenly increases dramatically. Meanwhile, the mean density profile in part of the domain inverts, suggesting the formation of a fully convective layer (see Section 5.2). This is reminiscent of the formation of gravity waves and layer development in oceanic fingering convection (Stellmach et al. 2011). We will discuss these large-scale instabilities and related transport properties in Section 5.

Zoom In Zoom Out Reset image size

To calculate the saturated fluxes in the homogenous phase of all simulations, i.e., prior to the onset of any large-scale dynamics discussed above, we use the following method. We identify the first local minimum in Nu_T(t) after saturation and set this to be the beginning of the saturated regime. We then fit the remainder of the time series to a linear function of time with a least-squares fit, which yields a slope and a y-intercept. We use regression analysis to estimate an uncertainty for each parameter. If the uncertainty in the slope is greater than its magnitude, we then say that the Nusselt curve after saturation is flat and set the end of the saturated regime at the end of the simulation. On the other hand, for a simulation like the one shown in Figure 3, where the flux saturates initially but gravity waves begin enhancing the transport again thereafter, we are only interested in the short interval of time just after saturation. If the linear fit is not sufficiently flat, we reduce the end time of the "homogeneous phase" progressively until the fitted slope is zero (within its uncertainty). Once the time interval for homogenous fingering convection has been identified, we calculate the saturated fluxes by averaging them over this interval. The error is taken to be the rms of the fluxes.

A compilation of our new results with published data from Traxler et al. (2011a) and Radko & Smith (2012) is given in Table 1. Figure 4 shows the corresponding values of the mean Nusselt numbers of both composition and temperature, in the homogenous saturated phase, as a function of the reduced density ratio r (see Equation (11)). As expected, we see that both Nu_T − 1 and Nu_μ − 1 tend to 0 as r → 1. In contrast, as r → 0, both Nu_T − 1 and Nu_μ − 1 increase rapidly, presumably approaching values of Rayleigh–Benard convection in a triply periodic domain (Calzavarini et al. 2005; Garaud et al. 2010). As an additional note, it is clear from Figure 4 that Nu_μ does not depend on Pr and τ individually, but rather on their ratio. This effect had already been observed in Traxler et al. (2011a), and led them to propose their empirical scaling laws:

$$\begin{eqnarray} \hbox{Nu}_T-1&=\hbox{Pr}^{1/2} \tau ^{3/2}f(r), \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} \hbox{Nu}_{\mu }-1&=\sqrt{\frac{\hbox{Pr}}{\tau }}g(r), \end{eqnarray} \tag{ 17 }$$

where g(r) and f(r) have the form ae^−br(1 − r)^c, where for f(r), a = 264 ± 1, b = 4.7 ± 0.2, c = 1.1 ± 0.1, and for g(r), a = 101 ± 1, b = 3.6 ± 0.3, c = 1.1 ± 0.1.

Zoom In Zoom Out Reset image size

As discussed in Section 1, the prescription given by Traxler et al. (2011a) agrees well with the data for high R₀; however, these scalings underestimate the transport at low R₀, particularly as Pr and τ decrease. It is clear from Figure 4 that their prescription breaks down at very low r, particularly as Pr decreases. For this reason we now propose a theory that better models the data at both large and small r. We derive this new theory with physical justification.

The properties of the fastest growing fingering mode can be determined by linearizing the governing equations and assuming exponential solutions of the form

$$\begin{equation} q=\hat{q}e^{\lambda {t}+i\mathbf {k}\cdot {\mathbf {x}}}, \end{equation} \tag{ 18 }$$

for each of the velocity, temperature, composition, and pressure perturbations (Baines & Gill 1969). This yields a cubic equation for the growth rate λ, with coefficients that depend on the wave vector, k, and the governing non-dimensional parameters, Pr, τ, and R₀. It can be shown that all properties of the fastest growing mode (e.g., the velocity, temperature, and salinity fields) are independent of z, hence the term "elevator" mode used to describe them. Since our equations are symmetric in x and y, λ only depends on the magnitude l of the horizontal wavenumber and not on its direction.

The resulting cubic then reads (Baines & Gill 1969):

$$\begin{eqnarray} \lambda ^3&+a_2\lambda ^2+a_1\lambda +a_0=0, \quad \hbox{ where} \nonumber\\ a_2&=l^2(1+\hbox{Pr}+\tau), \nonumber \\ a_1&=l^4(\tau \hbox{Pr}+\hbox{Pr}+\tau)+\hbox{Pr}\left(1-\frac{1}{R_0}\right), \nonumber \\ a_0&=l^6\tau \hbox{Pr}+l^2\hbox{Pr}\left(\tau -\frac{1}{R_0}\right). \end{eqnarray} \tag{ 19 }$$

To identify the fastest growing mode, we maximize λ with respect to the wavenumber l. This yields a quadratic for λ:

$$\begin{eqnarray} a_2\lambda ^2&+a_1\lambda +a_0=0, \quad \hbox{where} \nonumber\\ a_2&=1+\hbox{Pr}+\tau , \nonumber \\ a_1&=2l^2(\tau \hbox{Pr}+\tau +\hbox{Pr}), \nonumber \\ a_0&=3l^4\tau \hbox{Pr}+\hbox{Pr}\left(\tau -\frac{1}{R_0}\right). \end{eqnarray} \tag{ 20 }$$

Equations (19) and (20) can be solved numerically simultaneously to determine the wavelength and growth rate of the fastest growing mode. More interestingly for anyone interested in quick analytical estimates, it can be shown that in the limit of Pr, τ ≪ 1,

$$\begin{equation} \lambda \approx \left\lbrace \begin{array}{@{}ll}\sqrt{\hbox{Pr}} & \hbox{if r\ll \hbox{Pr}\ll 1} \\ \ms \sqrt{\frac{\hbox{Pr}\tau }{r}} & \hbox{if \hbox{Pr}\ll {r}\ll 1} \end{array} \right.\!\!, \end{equation} \tag{ 21 }$$

and

$$\begin{equation} l^2\approx \frac{1}{\sqrt{1+\tau /\hbox{Pr}}}. \end{equation} \tag{ 22 }$$

A detailed derivation of these asymptotic limits (and higher-order terms) is presented in Appendix B.

Radko & Smith (2012) consider all possible sources of secondary instabilities for the elevator modes. Here, for simplicity, we restrict our analysis to instabilities arising from the shear between adjacent elevators and neglect viscosity. As we demonstrate below, this yields very satisfactory results, and allows for a much simpler solution to the problem.

The elevator modes are characterized by a vertically invariant flow field of the kind

$$\begin{equation} \mathbf {U}(x,t)= w_0 e^{\lambda t} \sin (lx)\mathbf {e}_z \equiv w_{\rm E}(t) \sin (lx) \mathbf {e}_z\mbox{, } \end{equation} \tag{ 23 }$$

where λ and l are now strictly defined as the growth rate and horizontal wavenumber of the fastest growing mode, introduced in the previous section. Without loss of generality, we have selected the horizontal phase of the mode to be sin (lx). This flow field is associated with a temperature and composition field

$$\begin{equation} T_{\rm E}(x,t)=T_0 e^{\lambda t} \sin (lx) \quad \mbox{ and } \mu _{\rm E}(x,t)=\mu _0 e^{\lambda t} \sin (lx) \mbox{, } \end{equation} \tag{ 24 }$$

where w₀, T₀, and μ₀ are related via Equations (8) and (9):

$$\begin{eqnarray} &&\lambda T_0 + w_0 = - l^2 T_0 \Rightarrow T_0 = - \frac{w_0}{\lambda + l^2} \mbox{, } \end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray} &&\lambda \mu _0 + R_0^{-1} w_0 = - \tau l^2 \mu _0 \Rightarrow \mu _0 = - \frac{R_0^{-1} w_0}{ \lambda + \tau l^2} \mbox{, } \end{eqnarray} \tag{ 26 }$$

since elevator modes are exact solutions of the full set of nonlinear equations. At early times, the velocity shear, w_E(t), is weak and the elevator modes grow unimpeded. Using the definitions of the Nusselt numbers given in Equation (12) with Equations (23)–(25), we find that

$$\begin{eqnarray} \hbox{Nu}_T(t) &=& 1 - \frac{1}{L_xL_yL_z} \int _x \int _y \int _z w_{\rm E}(t) T_{\rm E}(t) \sin ^2(lx) dxdydz \nonumber \\ &=& 1 - \frac{1}{2} w_{\rm E}(t) T_{\rm E}(t) = 1 + \frac{w^2_0}{2(\lambda + l^2)} e^{2\lambda t} \mbox{.} \end{eqnarray} \tag{ 27 }$$

A similar expression can be obtained for Nu_μ(t). This shows that the Nusselt numbers initially grow exponentially with a growth rate that is twice that of the fastest growing mode, as mentioned in Section 3.1, prior to saturation.

We now consider secondary shearing instabilities on the elevator mode flow. Since Pr ≪ 1, we assume that viscosity is negligible. In that case, the growth rate of the shearing modes, σ, is a simple increasing function of the velocity within the fingers. Since the latter increases exponentially, σ also increases, while the growth rate of the fingering instability remains constant. We therefore expect that there will be a time where the two are of the same order of magnitude. At this point, the elevator modes are destroyed and saturation occurs.

As in Radko & Smith (2012), we neglect the temporal variation of the primary elevator modes when evaluating the growth rate of the secondary shearing modes, and simply write

$$\begin{equation} \mathbf {U}(x,z)= \hat{w}_{\rm E} \sin (lx)\mathbf {e}_z \mbox{.} \end{equation} \tag{ 28 }$$

The growth rate of the shearing instability can be found through dimensional analysis or through Floquet theory (see Appendix A). The dimensional argument goes as follows: the only relevant velocity in the system is that of fluid within the fingers, $\hat{w}_{\rm E}$, and the only relevant length scale is 1/l, so the only relevant growth rate is $\hat{w}_{\rm E}l$. Hence,

$$\begin{equation} \sigma = K \hat{w}_{\rm E}l \mbox{,} \end{equation} \tag{ 29 }$$

where K is a universal constant of proportionality. By "universal," we imply that this constant is independent of any system parameter (τ, Pr, R₀) or any property of the elevator mode. All information about the latter is contained in l and $\hat{w}_{\rm E}$.

Assuming that saturation occurs when the shearing instability growth rate is of the order of the fingering growth rate can be written mathematically as

$$\begin{equation} \sigma = K \hat{w}_{\rm E}l = c \lambda , \end{equation} \tag{ 30 }$$

where c is independent of the fundamental parameters of the system. This equation uniquely determines the velocity within the fingers at saturation to be

$$\begin{equation} \hat{w}_{\rm E} = \frac{C\lambda \sqrt{2}}{l}, \end{equation} \tag{ 31 }$$

where $C = c/K\sqrt{2}$ is another universal constant (to be determined). Equation (31) is similar to the one used by Denissenkov (2010) to estimate the effective diffusivity by dimensional analysis. Our own model, however, departs from his analysis at this point and produces a result that more accurately fits the results of simulations.

A given elevator mode with velocity $w_{\rm E} (x,z) = \hat{w}_{\rm E} \sin (lx)$ has a temperature profile $T_{\rm E} (x,z) = \hat{T}_{\rm E} \sin (lx)$ with $\hat{T}_{\rm E} = - ({\hat{w}_{\rm E}}/{\lambda + l^2})$, for the same reasons that led us to Equation (25). Hence, at saturation,

$$\begin{equation} {\rm Nu}_T = 1 + \frac{\hat{w}^2_{\rm E}}{2(\lambda + l^2)} = 1 + C^2 \frac{ \lambda ^2 }{l^2 (\lambda + l^2)} \end{equation} \tag{ 32 }$$

using Equation (31). And by similar arguments, it can be shown that

$$\begin{equation} \hbox{Nu}_{\mu }=1 + C^2\frac{\lambda ^2}{\tau l^2(\lambda +\tau {l^2})}, \end{equation} \tag{ 33 }$$

where C is the same constant as in Equation (32). The physical interpretation of C is now clearer. Since C is related to the relative importance of the shearing and fingering instabilities at saturation, a large value of C indicates that the shearing instability cannot easily disturb the initial fingers, resulting in a large saturated flux. Conversely, a small value suggests that shear easily destroys fingers, resulting in a small saturated flux.

We compare this prescription with data from this study and those simulations from Traxler et al. (2011a) and Radko & Smith (2012) that approach the astrophysical regime. We fit C by using a chi-squared statistical test and find a best fit at C = 7.0096. As illustrated in Figure 5, for a value of C = 7 and τ ⩽ Pr, we find that our theoretical results match all simulations remarkably well. It is important to note that, while C needs to be fitted, it is a universal constant and cannot depend on Pr, τ, or R₀. The fact that we are able to use only one value of C to correctly fit the data suggests that the premises of this theory are correct.

Zoom In Zoom Out Reset image size

In Figure 6, we show that in the opposite case (τ > Pr) we do not observe the same quality of fit. The full analysis by Radko & Smith (2012) is also unable to explain these results, so it is likely that saturation in this case operates somewhat differently. However, since stars are always in the former situation, the poorness of fit of the latter case is interesting but not relevant to our ultimate purpose.

Zoom In Zoom Out Reset image size

As such, Equations (32) and (33) yield the thermal and compositional Nusselt numbers provided λ and l are known. Calculating these quantities requires the simultaneous solution of a cubic and a quadratic (see Section 4.1), which can be done numerically quite easily. However, analytical approximations to Nu_T and Nu_μ can also be obtained using asymptotic analysis (see Appendix B), to find that for r ≪ (τ, Pr) ≪ 1,

$$\begin{eqnarray} \hbox{Nu}_T-1&\approx {C}^2\hbox{Pr}\left(1+\frac{\tau }{\hbox{Pr}}\right), \end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray} \hbox{Nu}_{\mu }-1&\approx {C}^{2}\sqrt{\frac{\hbox{Pr}}{\tau }}\sqrt{1+\frac{\tau }{\hbox{Pr}}}. \end{eqnarray} \tag{ 35 }$$

For the case with (τ, Pr) ≪ r ≪ 1 and τ/Pr ∼ 1,

$$\begin{eqnarray} \hbox{Nu}_T-1&\approx {C^2}\hbox{Pr}\tau \frac{1+\frac{\tau }{\hbox{Pr}}}{r}, \end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray} \hbox{Nu}_{\mu }-1&\approx {C^2}\sqrt{\frac{\hbox{Pr}}{\tau }}\sqrt{\frac{1+\frac{\tau }{\hbox{Pr}}}{r}}. \end{eqnarray} \tag{ 37 }$$

Note that in this case, Nu_μ depends only on τ/Pr and not on Pr or τ individually, as was noticed by Traxler et al. (2011a). And finally for the case as r → 1, (τ, Pr) ≪ 1,

$$\begin{eqnarray} &&\hbox{Nu}_T-1\approx {C^{2}}\tau ^{2}\frac{4}{9}\left(1-r\right)^{2}, \end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray} &&\hbox{Nu}_{\mu }-1\approx {C^{2}}\frac{4}{9}\left(1-r\right)^{2}. \end{eqnarray} \tag{ 39 }$$

Let us discuss these results in more detail. First note that in the case (τ, Pr) ≪ r ≪ 1, we find Nu_T − 1∝Prτ and Nu_μ − 1∝Pr^1/2τ^−1/2, and for r → 1, we find Nu_T − 1∝τ² while Nu_μ − 1 has no pre-factor dependence on Pr or τ. These scalings are reasonably (though not exactly) consistent with the results of Traxler et al. (2011a), who found Nu_T − 1∝Pr^1/2τ^3/2 and Nu_μ − 1∝Pr^1/2τ^−1/2. The discrepancy is not so surprising given that they did not differentiate between different regimes. Second, note that in astrophysical cases, Pr ≪ 1, which implies that turbulent heat transport by homogeneous fingering convection is negligible. Such is not the case for Nu_μ − 1, which can become quite large, especially in the limit of r → 0. We will discuss the implications of these results more in Section 6. Finally, note that these asymptotic formulae should work well for the true astrophysical regime where Pr, τ ≪ 1. However, since none of our simulations are particularly close to this regime, the asymptotic approximation does not compare well with the numerical simulations we have run so far.

Radko (2003) developed the so-called γ-instability theory to explain the formation of thermohaline staircases in fingering regions of the ocean. The theory was found to agree with direct numerical simulations of staircase formation in two and three dimensions (Radko 2003; Stellmach et al. 2011) and has since become well accepted in the physical oceanography community.

The original γ-instability theory is a linear mean-field instability theory. Radko (2003) first averaged the governing equations of fingering convection (see Equations (6)–(9)) spatially, over length scales that span several fingers, and temporally, over multiple finger lifetimes. He then performed a linear stability analysis of the averaged, or "mean-field," equations. He argued that the stability of the system depends on the quantity γ_turb, defined as the ratio of the turbulent thermal flux, Nu_T − 1, to the turbulent compositional flux, τ(Nu_μ − 1)/R₀. More specifically, he showed that homogeneous fingering convection is unstable to layering if and only if ∂γ_turb/∂R₀ < 0. Traxler et al. (2011a) later improved this theory and showed that in the astrophysical case, this result still holds as long as γ is redefined as the ratio of the total fluxes, diffusive and turbulent⁷; in other words, a system is unstable to layering if and only if

$$\begin{equation} \frac{\partial {\gamma }}{\partial {R_{0}}} < 0 \quad \hbox{ with } \gamma =\frac{R_{0}\hbox{Nu}_T}{\tau \hbox{Nu}_{\mu }}. \end{equation} \tag{ 40 }$$

As mentioned in Section 3.2, some of our simulations at very low r exhibit properties consistent with layering. In order to determine whether these layers form via the γ-instability, we plot in Figure 7 the variation of γ with R₀ for various values of Pr and τ. Unlike γ_turb, which does decrease for some r at low Pr, τ, the quantity γ always seems to increase. As discussed by Traxler et al. (2011a), this is because Nu_T and Nu_μ are both dominated by diffusive fluxes at low Pr and τ. The diffusive flux ratio, R₀/τ, increases strongly with R₀ because τ is asymptotically small and R₀ approaches 1/τ as the system becomes increasingly stable. This effect dominates the total flux ratio. Since γ is always an increasing function of R₀, we conclude that the observed layers cannot be forming by the γ-instability. Thus, while the γ-instability can lead to the development of staircases in high Pr oceanic fingering, it is unlikely to be relevant for astrophysical fingering, confirming the results of Traxler et al. (2011a).

Zoom In Zoom Out Reset image size

Despite the fact that γ is a strictly increasing function of r, we do observe layer formation in some simulations. We now study these results in more detail, focusing on the case presented in Figure 3, which has Pr = 1/10, τ = 1/30, and R₀ = 1.1 and shows strong evidence for the formation of a convective layer. Figure 8 shows a snapshot of the compositional perturbation at t = 825; large-scale plumes characteristic of overturning convection are clearly visible. We can check to see whether these plumes are associated with an inversion of the horizontally averaged total density profile, which is given by

$$\begin{equation} \langle \rho _{{\rm tot}}\rangle =-\left(1-R_0^{-1}\right)z-\langle {T}\rangle +\langle {\mu }\rangle. \end{equation} \tag{ 41 }$$

Figure 9 shows 〈ρ_tot〉 as a function of height for the simulation shown in Figure 3 at four different times selected during the later stages, where the Nusselt numbers increase significantly. The solid line indicates the linear background gradient for reference. Note how each of the four density profiles shows a sharp, stably stratified transition region between z = 25 and z = 45, the interface. Three of the four profiles also have a region where the total density increases with z, the convective layer. This demonstrates that the homogeneous fingering has indeed transitioned into a layered system with convective layers separated by stable interfaces. At t = 725, shown in Figure 9, the convective layer is just beginning to develop, and a slight inversion of the density profile appears in Figure 9. As the simulation approaches the stage of maximum transport at t = 750 (see Figure 3), the density inverts more strongly. The ensuing period of active mixing thickens the interface again and flattens the density profile. At t = 775, the convective mixing briefly shuts down. This process repeats, as can be seen at t = 800, when the density profile once more is beginning to invert. Note that the interface appears to move up and down during this period, which we attribute to the advection of a large quantity of compositionally dense material by convective plumes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This is not the only such simulation in which we observe the formation of staircases: we see this behavior at low r (≲ 0.003) for several values of Pr and τ, including Pr = 1/10, τ = 1/10 and Pr = 1/3, τ = 1/10. In all these cases, the Nu_T and Nu_μ time series as well as the properties of the layered state are qualitatively similar to the ones shown in Figures 3, 8, and 9.

Having established that this layering transition cannot be due to the γ-instability, we look instead at the original mechanism proposed by Stern (1969) (see also Stern et al. 2001) for the formation of thermohaline staircases in the ocean. In addition to the γ-instability, fingering convection is also known to excite large-scale gravity waves through another mean-field instability, the "collective instability." Stern (1969) studied this effect, and showed that large-scale gravity waves can be excited under some circumstances and grow in amplitude exponentially. These waves transport momentum, heat, and heavy elements until they begin to "break," at which point the advected material mixes with the background. He went on to suggest that if the waves break on large enough scales, this could cause enough local mixing to trigger the formation of staircases.

In related numerical simulations at oceanographically relevant parameters (Pr = 7, τ = 0.3), Stellmach et al. (2011) found that gravity waves are indeed excited by the collective instability but do not trigger layers. Layer formation in their simulations—and thus presumably in the ocean—is caused by the γ-instability rather than the collective instability. Our simulations suggest, however, that the converse may be true in the astrophysical parameter regime. As discussed in Figure 3, gravity waves are observed to dominate the dynamics of the system between t = 100 and t = 600. Furthermore, the gradual increase in Nu_T and Nu_μ during that phase suggests that the waves have high enough amplitudes to interact nonlinearly with each other, and with the background, by contrast with the oceanographic case. Indeed, linear waves cause oscillations in Nu_T and Nu_μ without changing their mean values—only nonlinear waves can affect the mean. Since the waves in Figure 3 are clearly strong enough to interact nonlinearly, we can hope that they may also break and cause mixing on larger scales, following the original hypothesis proposed by Stern (1969).

Checking this hypothesis in detail is difficult with the available simulations. However, we can at least determine whether our observed waves do indeed appear in accordance with the collective-instability theory. A system is unstable to the collective instability if the Stern number, in our notation defined as

$$\begin{equation} A=\frac{(\hbox{Nu}_T-1) \left(\gamma _{{\rm turb}}^{-1}-1 \right)}{\hbox{Pr}\left(1-R_{0}^{-1}\right)}, \end{equation} \tag{ 42 }$$

exceeds order unity (Stern et al. 2001). Note that this expression uses γ_turb = R₀(Nu_T − 1)/τ(Nu_μ − 1) instead of γ. For the simulation shown in Figures 3 and 9, Nu_T − 1 = 6.55 and γ_turb = 0.438, so that the Stern number is 924. This implies that the system is indeed unstable to the collective instability.

More generally, we find that all the simulations in which we observe layer formation have Stern numbers of the order of 10³ or more. However, not all such simulations go into layers, as can be seen in Figure 10, so it is likely that the condition for layer formation through the collective instability also depends on the value of r. It is also possible that these gravity waves grow in all simulations with A > 1, but since the growth rate depends on the Stern number, gravity waves may be too weak to detect if A is too low. If this is the case, it will be difficult to see layer formation in models with smaller A because these simulations are expensive and difficult to integrate for long. We discuss this more in Section 6.3.

Zoom In Zoom Out Reset image size

Much work still needs to be done to identify the parameters for which we expect layer formation in the astrophysical regime and to quantify transport in the layered case (see Wood et al. 2012 for related efforts in the case of semi-convection). We have only three simulations that develop into layers and thus require a more extensive sample of simulations at low r in order to better understand the conditions of layer formation and their properties as a function of Pr, τ, and r. In every simulation that exhibits layers, we find that a single convective plume spans the majority of the domain and do not see the development of multiple layers, as normally seen in simulations of the fingering regime in the oceanographic case (e.g., Stellmach et al. 2011) and in the case of semi-convection (Rosenblum et al. 2011; Mirouh et al. 2012). To characterize this process, we must increase the size of the domain and extend the integration time to let the layers develop more naturally. This will allow us to describe more completely the size scales and transport of thermocompositional layers in astrophysical fingering convection. This problem, however, is beyond the scope of this paper and is deferred to a later publication.

The first step toward using our model in stellar evolution calculations is to identify the non-dimensional parameters of the system in question at each grid point and/or time step. Recall that Pr = ν/κ_T and τ = κ_μ/κ_T. The expressions for the microscopic viscosity and diffusivities can be found in many textbooks of plasma physics (e.g., Chapman & Cowling 1970). In stars, the value of Pr typically ranges from 10⁻⁷ to 10⁻⁶ and that of τ from 10⁻⁸ down to 10⁻⁶, with Pr > τ.

The density ratio, R₀, is given by

$$\begin{equation} R_0=\frac{\alpha {\left(T_{0z}-T_{0z}^{{\rm ad}}\right)}}{\beta {\mu _{0z}}}=\frac{\nabla -\nabla _{{\rm ad}}}{\frac{\phi }{\delta }\nabla _\mu }, \end{equation} \tag{ 43 }$$

where ∇ = (dln T/dln P) is the actual temperature gradient, ∇_ad = (dln T/dln P)_S is the adiabatic temperature gradient, ∇_μ = (dln μ/dln P) is the mean molecular weight gradient, and δ = −(∂ln ρ/∂ln T)_{P, μ} and ϕ = (∂ln ρ/∂ln μ)_{P, T} (see Kippenhahn & Weigert 1990 for a detailed description of these gradients). The quantities ∇ and ∇_μ can be calculated from the stellar model grid. The other values, ∇_ad, ϕ, and δ, can be calculated from the equation of state. Note that R₀ can take a broad range of values in stellar interiors, but the fluid is fingering-unstable only when 1 < R₀ < τ⁻¹.

Once Pr, τ, and R₀ are known, one can then calculate Nu_T and Nu_μ according to Equations (32) and (33) for a numerically precise answer or by the asymptotic expressions in Section 4.3 if a rough estimate is all that is needed. Since Nu_μ is the total flux of composition divided by the background diffusive flux, the effective compositional diffusivity for modeling transport by fingering convection is given by

$$\begin{equation} D_{\mu }=\hbox{Nu}_{\mu }\kappa _{\mu }, \end{equation} \tag{ 44 }$$

where

$$\begin{equation} F_\mu =-D_{\mu }\mu _{0z} \end{equation} \tag{ 45 }$$

is the total dimensional flux of the mean molecular weight through the fingering region. An analogous expression can be derived for the heat transport, but this term is nearly always negligible (see Section 4.2) when Pr ≪ 1 unless layers form. We ignore it from here onward. By contrast, the turbulent compositional transport can be significant, as can be seen in Figure 11, where we plot the logarithm of Nu_μ − 1 as a function of τ/Pr = κ_μ/ν and r, assuming Pr = 10⁻⁶. Near the marginal stability limit, r = 1, the compositional turbulent transport is negligible as well, and Nu_μ − 1 drops to zero. However, for low r, and particularly for low τ/Pr, the turbulent compositional transport increases by many orders of magnitude. This plot is intended as a quick estimate for a large range of parameters; more accurate estimates can be obtained from the expressions given in Sections 4.2 and 4.3.

Zoom In Zoom Out Reset image size

The results presented so far concern only the case of homogeneous fingering convection. As discussed in Section 5.2, however, layer formation seems to be possible in some cases and leads to a significant increase in transport compared with the homogeneous levels. To determine when this may occur, we use the theory from Section 4 to estimate the Stern number as a function of τ/Pr = κ_μ/ν and r and show the results in Figure 12. We find that the Stern number is largest for small r and large τ/Pr. Since we observe layer formation for a Stern number near or above 10³ (see Section 5.2), the same criterion suggests naively that layer formation will only be relevant in stellar interiors (where τ/Pr < 0.1) for r ≲ 10⁻⁵. This is so close to actual overturning convection that it may in fact be irrelevant for practical purposes. Of course, as mentioned in Section 5.2, much work still needs to be done to confirm these results.

Zoom In Zoom Out Reset image size

We now reexamine several problems that have been recently discussed in the context of fingering convection, namely, planet pollution (Vauclair 2004; Garaud 2011) and the impact of ³He fusion in red giant branch stars (Charbonnel & Zahn 2007; Denissenkov 2010).

As first discussed by Vauclair (2004), fingering convection determines the timescale during which planetary material remains detectable on the surface of a star after infall. Garaud (2011) used the fingering transport prescription from Traxler et al. (2011a) to study the problem and determined that this timescale is too short for the excess metallicity of planet-host stars to be related to planet infall. Since we have found that the prescription by Traxler et al. (2011a) only underestimates the transport at low r while still recovering the scalings at high r, the use of our new, improved model can only increase the turbulent mixing rates and decrease the mixing time, so the qualitative conclusion of Garaud (2011) still holds: the fact that planets are more readily found around high metallicity host stars must be a primordial effect.

To address the case of red giant branch stars, we compare our results to the numerical simulations of Denissenkov (2010). For his choice of governing parameters, Pr = 4 × 10⁻⁶, τ = 2 × 10⁻⁶, and R₀ = 1700, he finds Nu_μ = 998. He concludes from this result that fingering convection cannot provide all the additional mixing needed in red giant stars. Our model finds a remarkably similar value of Nu_μ = 1294 for the same parameters, and thus we confirm his conclusion that additional mixing is necessary to explain observations.

We have developed a simple semi-analytical model for the vertical transport of heat and composition by homogeneous fingering convection that reproduces the results of numerical simulations from this study and from previous studies and can easily be implemented in a stellar evolution code. We have found that, under conditions relevant for stellar interiors, the turbulent heat transport is negligible, but compositional transport can be important. In particular, this model predicts more efficient transport for weakly stratified systems than was initially suggested in previous work by Traxler et al. (2011a).

We have also found the first evidence for layer formation by fingering convection in the astrophysical parameter regime. We have identified the collective instability as the origin of layer formation. By contrast, in the oceanographic and semi-convective cases, the γ-instability causes layer formation (e.g., Stellmach et al. 2011; Mirouh et al. 2012). Much work still needs to be done to characterize the conditions for layer formation and transport by layered convection in this case. With our limited computational resources, we have found that layer formation appears to require high Stern numbers (see Equation (42)) to occur. If true, this would imply that layers are only likely to form in regions that are very close to being fully convective anyway. However, it remains to be seen whether lower Stern numbers can also result in layer formation. To study this will require additional simulations and longer integration times than have been covered in this paper. In addition to the conditions of layer formation, it is also important to develop a theoretical or empirical model for the transport of this case. Our simulations show that transport increases significantly when layers form. Concurrent work by Wood et al. (2012) reveals that thermal and compositional transport in the layered semi-convective case scales with the layer height to the power of 4/3. If this scaling also applies here, we may expect that very significant transport rates could occur for large enough layer heights. What determines the latter, however, will also require new theories.

In the limit of r ≪ (Pr, τ), we find numerically that λ is independent of r. The same is true for the wavenumber l of the fastest growing mode. This suggests that both λ and l are continuous in the limit r → 0, and that one may study this first regime simply by solving Equations (19) and (20) with r ≡ 0 (alternatively, R₀ = 1). The resulting equations are (using the definition of ϕ given above)

$$\begin{eqnarray} &&\lambda ^3 + \lambda ^2 l^2 (1 + \hbox{Pr}(1+\phi)) + \lambda l^4 \hbox{Pr}(1 + \phi (\hbox{Pr}+1))\nonumber\\ &&\quad + l^6 \hbox{Pr}^2 \phi + l^2 \hbox{Pr}(\hbox{Pr}\phi - 1) =0 \mbox{,} \end{eqnarray} \tag{ B1 }$$

and

$$\begin{eqnarray} &&(1 + \hbox{Pr}(1+\phi)) \lambda ^2 + 2 \lambda l^2 \hbox{Pr}(1 + \phi (\hbox{Pr}+1))\nonumber\\ &&\quad + 3 l^4 \hbox{Pr}^2 \phi + \hbox{Pr}(\hbox{Pr}\phi - 1) =0 \mbox{.} \end{eqnarray} \tag{ B2 }$$

Subtracting l² times the second equation from the first yields the simplified cubic:

$$\begin{equation} \lambda ^3 - \lambda l^4 \hbox{Pr}(1 + \phi (\hbox{Pr}+1)) - 2 l^6 \hbox{Pr}^2 \phi =0 \mbox{.} \end{equation} \tag{ B3 }$$

Next, we study the behavior of λ and l with varying Pr and τ (via ϕ). We observe that the exact solution for λ varies as Pr^1/2, while l² is more-or-less independent of Pr. Using this information, we propose the following asymptotic expansions:

$$\begin{eqnarray} &&\lambda = \hbox{Pr}^{1/2}(\lambda _0 + \lambda _1 \hbox{Pr}^{1/2}) + \cdots \quad \mbox{ and } l^2 = l_0^2 + l^2_1 \hbox{Pr}^{1/2} + \cdots\nonumber\\ \end{eqnarray} \tag{ B4 }$$

expanding l² rather than l since Equations (19) and (20) only depend on the former.

Plugging these two ansatz into Equations (B2) and (B3), and solving for all unknowns at their respective orders in Pr yields

$$\begin{eqnarray} & \lambda = \sqrt{\hbox{Pr}}-\hbox{Pr}\sqrt{1+\phi } + \cdots, \quad \mbox{and} \nonumber \\ & l^2 = \frac{1}{\sqrt{1+\phi }}-\sqrt{\hbox{Pr}}\left(1+\frac{\phi }{(1+\phi)^2}\right)+\cdots. \end{eqnarray} \tag{ B5 }$$

As can be seen in Figure 13, this approximation does well for small Pr and τ for r < (Pr, τ). Note that although the only illustrated cases are those with Pr = τ, we have tested cases down to ϕ ∼ 10⁻².

We may then use Equations (32) and (33) to estimate the thermal and compositional Nusselt numbers. Keeping only the leading order terms yields

$$\begin{eqnarray} \hbox{Nu}_T-1&\approx {C}^2\left(\hbox{Pr}\left(1+\phi \right)-\hbox{Pr}^{3/2}\frac{1+\phi ^2}{\sqrt{1+\phi }}+\cdots \right), \end{eqnarray} \tag{ B6 }$$

$$\begin{eqnarray} \hbox{Nu}_{\mu }-1&\approx {C^{2}}\left(\frac{\sqrt{1+\phi }}{\phi \sqrt{\hbox{Pr}}}-\frac{\phi }{1+\phi }+\cdots \right). \end{eqnarray} \tag{ B7 }$$

This too is compared to the numerical solution in Figure 14.

Zoom In Zoom Out Reset image size

In this regime, since (Pr, τ) ≪ r, we expand Equations (19) and (20) first in Pr and then in r. Inspection of the numerical solutions suggests that λ is proportional to Pr, while l² is more-or-less independent of it. Seeking solutions of Equations (19) and (20) of the form $\lambda = \hbox{Pr}\hat{\lambda }$, to the lowest order in Pr, yields

$$\begin{eqnarray} l^2 \hat{\lambda }^2+ [l^4\left(1+\phi \right)+1] \hat{\lambda }+l^6\phi +l^2\phi \frac{r-1}{r}&=0, \quad \mbox{ and}\quad \end{eqnarray} \tag{ B8 }$$

$$\begin{eqnarray} \hat{\lambda }^2+2l^2\left(1+\phi \right)\hat{\lambda }+3l^4\phi +\phi \frac{r-1}{r}&=0. \end{eqnarray} \tag{ B9 }$$

Eliminating the $\hat{\lambda }^2$ term yields

$$\begin{equation} \hat{\lambda }[ l^4(1+\phi)-1]+2l^6\phi =0 \mbox{,} \end{equation} \tag{ B10 }$$

which can be easily solved to find

$$\begin{equation} \hat{\lambda }= -\frac{2l^6\phi }{l^4(1+\phi)-1} \mbox{.} \end{equation} \tag{ B11 }$$

Plugging this into Equation (B9) and collecting the results in terms of l yields

$$\begin{eqnarray} &&- l^{12} (1-\phi)^2 + l^8 (1+\phi) \left[ \phi - 1 - \frac{1+\phi }{r} \right]\nonumber\\ &&\quad + l^4 \left[ 1 - 2\phi + \frac{2}{r} (1+\phi) \right] + \frac{r-1}{r} = 0, \end{eqnarray} \tag{ B12 }$$

which is a cubic equation for l⁴. Unfortunately, this cubic has no simple solution, so we look for an expansion in the other asymptotic variable, r. Inspection of the numerical solutions suggests that the expansion needs to be in terms of r^1/2, so we set $l^4=l_0^4\,+\,l_1^4\sqrt{r}\,+\,l_2^4 r + \cdots$. Using this ansatz in Equation (B12) yields, to lowest order,

$$\begin{equation} l_0^8(1+\phi)^2 -2 l_0^4 (1+\phi) + 1 = 0, \end{equation} \tag{ B13 }$$

which we can solve to find $l_0^4=1/(1+\phi)$. This, on its own, causes $\hat{\lambda }$ to diverge (see Equation (B11)), so we need to go to the next order in $\sqrt{r}$. By choosing the growing solution, we then find that

$$\begin{equation} l^4=\frac{1}{1+\phi } - 2\frac{\sqrt{r\phi }}{(1+\phi)^{5/2}}+\cdots , \end{equation} \tag{ B14 }$$

which eventually yields

$$\begin{equation} \lambda =\hbox{Pr}\sqrt{\frac{\phi }{r}}-\hbox{Pr}\sqrt{1+\phi }+\cdots , \end{equation} \tag{ B15 }$$

keeping the first two terms in the expansion only. As expected from the numerical solutions discussed above, λ scales like r^−1/2Pr to lowest order.

As in the previous section, we use Equations (32) and (33) along with these expansions to find approximate expressions for Nu_T and Nu_μ:

$$\begin{eqnarray} \hbox{Nu}_T-1&\approx {C^2}\hbox{Pr}^2\left(\phi \frac{1+\phi }{r}-2\sqrt{\frac{\phi +\phi ^2+\phi ^3}{r\left(1+\phi \right)}}+\cdots \right), \quad \mbox{and}\nonumber\\ \end{eqnarray} \tag{ B16 }$$

$$\begin{eqnarray} \hbox{Nu}_{\mu }-1&\approx {C^2}\left(\sqrt{\frac{1+\phi }{r\phi }}-\frac{1+2\phi +2\phi ^2}{\phi \left(1+\phi \right)}+\cdots \right). \qquad \end{eqnarray} \tag{ B17 }$$

We compare the results of the asymptotic expansion to the numerical solution in Figures 13 and 14, where we look at λ and Nu_μ − 1, respectively. We choose to set the transition between the two regimes at r = (τ, Pr). For those simulations where τ is several orders of magnitude smaller than Pr, we recommend choosing this point at r = τ. It is clear that the asymptotic expansion does an excellent job of reproducing the true numerical solution in the limit considered. We also see that, as suggested by the numerical data, Nu_μ only depends on ϕ rather than on Pr or τ individually as was noticed by Traxler et al. (2011a).

Equations (B11) and (B12) were derived without any assumptions on the value of r, aside from the fact that it needs to be much larger than Pr and τ. We can therefore use them directly to study the limit of r → 1. To do so, we set = 1 − r, and rewrite Equation (B12) in terms of the small parameter :

$$\begin{eqnarray} &&l^{12} (1-\phi)^2\left(\epsilon -1\right) - l^8(1+\phi) \left[ 2 + \epsilon (\phi -1)) \right]\nonumber\\ &&\quad + l^4 \left[3+\epsilon \left(2\phi -1\right)\right] - \epsilon = 0. \end{eqnarray} \tag{ B18 }$$

Inspection of the solutions for l⁴ in the limit of small suggests the expansion $l^4 = \epsilon (l_0^4 +\epsilon l_1^4 + \cdots)$. Plugging this into Equation (B18) then yields, to lowest order,

$$\begin{equation} l^4 = \frac{1}{3} (1-r)+\cdots. \end{equation} \tag{ B19 }$$

Plugging this expression into Equation (B11), we get

$$\begin{eqnarray} &&\lambda = 2 \hbox{Pr}\phi \left[ \frac{1}{3} (1-r) +\cdots \right]^{3/2}\left[ 1- \frac{1+\phi }{3} (1-r) +\cdots \right]^{-1}.\nonumber\\ \end{eqnarray} \tag{ B20 }$$

As before, we use Equations (B19) and (B20) with Equations (32) and (33) to find expansions in Pr and r for Nu_T and Nu_μ in the limits Pr → 0 and r → 1. We find

$$\begin{eqnarray} \hbox{Nu}_T-1&\approx {C^{2}}\hbox{Pr}^{2}\!\left(\!\frac{4}{9}\phi ^{2}\left(1-r\right)^{2}-\frac{8}{81}\phi ^{2}\left(\phi -8\right)\left(1-r\right)^{3}\cdots\!\! \right), \nonumber \\ \hbox{Nu}_{\mu }-1&\approx {C^{2}}\left(\frac{4}{9}\left(1-r\right)^{2}-\frac{8}{81}\left(\phi -5\right)\left(1-r\right)^{3}+\cdots \right).\nonumber \\ \end{eqnarray} \tag{ B21 }$$

We compare the asymptotic expressions to the numerical solutions for λ and Nu_μ in Figures 15 and 16. The expansion clearly works best for r ≳ 0.5. Note that these asymptotic expansions are designed to work only when ϕ ∼ 1 and should be used with caution if ϕ ≪ r.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size