Mixing via Thermocompositional Convection in Hybrid C/O/Ne White Dwarfs

We focus on a small region located in the vicinity of the CO/ONe core–mantle interface. We use a 2D Cartesian coordinate system (x, z), with gravity given by ${\boldsymbol{g}}=-g{{\boldsymbol{e}}}_{z}$. We are restricted to 2D simulations because of the vast range of scales that need to be accounted for numerically—from the tiniest scale associated with double-diffusive fingering, to the global scale associated with convection. This assumption is discussed in more detail in Section 5. We assume for simplicity that the temperature has a linear background state and the mean molecular weight profile has a constant background state

$$\begin{eqnarray}&&{T}_{0}(z,t)={T}_{{\rm{m}}}+{{zT}}_{0z}(t)\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\mu }_{0}(z,t)={\mu }_{{\rm{m}}}.\end{eqnarray} \tag{ 5 }$$

The values T_m and μ_m given here can be interpreted as the mean temperature and mean molecular weight per free electron of the pure CO core. The unstable compositional profile due to the CO/ONe transition, as well as perturbations to that state, will be added as initial conditions to drive the thermocompositional instability (see Section 3.2). While T_m could be time-dependent in a real WD, that time-dependence does not affect the hydrodynamics of the system much, so we neglect it here as a first approximation. By contrast, the temperature gradient T_0z is assumed to be time-dependent to model the effect of cooling, as in

$$\begin{eqnarray}&&{T}_{0z}(t)={T}_{0z}(t=0)\exp (-t/{\tau }_{\mathrm{cool}}),\end{eqnarray} \tag{ 6 }$$

where τ_cool is a characteristic global cooling timescale. This is a simple way of accounting for the gradual flattening of the temperature profile toward an isotherm, and the associated decay of the stabilizing temperature stratification, as the star cools.

The total pressure in the WD material can be approximated as a zero-temperature electron gas plus an ideal gas of ions

$$\begin{eqnarray}&&P={P}_{{\rm{e}}}+{P}_{\mathrm{ion}}=K{\left(\displaystyle \frac{\rho }{{\mu }_{{\rm{e}}}}\right)}^{\gamma }+\displaystyle \frac{\rho {RT}}{\bar{A}},\end{eqnarray} \tag{ 7 }$$

where we recall that ${\mu }_{{\rm{e}}}$, the mean molecular weight per free electron, is defined as ${\mu }_{{\rm{e}}}\equiv 1/{Y}_{{\rm{e}}}=\bar{A}/\bar{Z}$. This assumes a simple polytropic approximation for the electron pressure, with a constant K and an adiabatic index γ, which would fall in the range 5/3 (nonrelativistic electrons) to 4/3 (relativistic electrons). At the conditions of interest, the ion pressure contributes at the percent level. Finite-temperature corrections to the electrons and nonideal contributions from the ions, both of which are neglected in the above expression, also enter at the percent level. In what follows, we use the Boussinesq approximation for gases (Spiegel & Veronis 1960), in which fluid parcels are assumed to remain in pressure equilibrium with their surroundings at all times. Linearizing the equation of state (Equation (7)), and setting pressure perturbations δP to zero, we obtain an explicit expression for density perturbations δρ in terms of temperature perturbations δT, and perturbations to the mean molecular weight per free electron δμ_e:

$$\begin{eqnarray}&&\displaystyle \frac{\delta \rho }{\rho }\equiv -\alpha \delta T+\beta \delta {\mu }_{{\rm{e}}}\approx -\left(\displaystyle \frac{{P}_{\mathrm{ion}}}{\gamma P}\right)\displaystyle \frac{\delta T}{T}+\displaystyle \frac{\delta {\mu }_{{\rm{e}}}}{{\mu }_{{\rm{e}}}},\end{eqnarray} \tag{ 8 }$$

where we have used the fact the ion pressure is a small fraction of the total in order to neglect its composition dependence. This expression looks like the expression for an ideal gas, but with a thermal expansion coefficient, α, that is suppressed by a factor of P_ion/(γ P) ∼ 0.01, and the mean molecular weight per free electron ${\mu }_{{\rm{e}}}$ playing the analogous role to the total mean molecular weight. Recognizing this analogy, for notational convenience we will drop the subscript "e" and elide the "per free electron" in what follows.

Then, the equations for the perturbations around the mean state ($T={T}_{0}+\widetilde{T}$, $\mu ={\mu }_{0}+\widetilde{\mu }$) satisfy

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{u}}}{\partial t}+{\boldsymbol{u}}\cdot {\rm{\nabla }}{\boldsymbol{u}}=-\displaystyle \frac{{\rm{\nabla }}p}{{\rho }_{{\rm{m}}}}+g(\alpha \widetilde{T}-\beta \widetilde{\mu }){{\boldsymbol{e}}}_{z}+\nu {{\rm{\nabla }}}^{2}{\boldsymbol{u}},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \widetilde{T}}{\partial t}+{\boldsymbol{u}}\cdot {\rm{\nabla }}\widetilde{T}+w\left({T}_{0z}(t)-{T}_{\mathrm{ad},z}\right)={\kappa }_{T}{{\rm{\nabla }}}^{2}\widetilde{T},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \widetilde{\mu }}{\partial t}+{\boldsymbol{u}}\cdot {\rm{\nabla }}\widetilde{\mu }={\kappa }_{\mu }{{\rm{\nabla }}}^{2}\widetilde{\mu },\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{u}}=0\end{eqnarray} \tag{ 12 }$$

where ${\boldsymbol{u}}=(u,v,w)$ is the fluid velocity, p is the pressure, ρ_m is the mean density of the region, ν is the kinematic viscosity, κ_T is the thermal diffusivity, κ_μ is the compositional diffusivity, and ${T}_{\mathrm{ad},z}=-g/{c}_{P}$ is the background adiabatic temperature gradient. The relevant compositional diffusivity is the diffusivity of the ions, as charge neutrality means that the mean molecular weight per free electron can only change due to the transport of ions. The relevant thermal diffusivity is dominated by electron transport. The thermal expansion and compositional contraction coefficients are denoted as α and β, respectively. All of these quantities are assumed to be constant for simplicity. This assumption can of course affect detailed predictions of the model, but not its general conclusions (see the discussion in Section 5).

We then nondimensionalize the equations. We take as the characteristic lengthscale $[l]=d={({\kappa }_{T}\nu /(g\alpha | {T}_{\mathrm{ad},z}| ))}^{1/4}$. The characteristic time, temperature, and composition scales are $[t]={d}^{2}/{\kappa }_{T}$, $[T]=d| {T}_{\mathrm{ad},z}|$, and $[\mu ]=(\alpha /\beta )d| {T}_{\mathrm{ad},z}|$, respectively. We define the Prandtl number to be $\Pr \equiv \nu /{\kappa }_{T}$, and the diffusivity ratio to be $\tau \equiv {\kappa }_{\mu }/{\kappa }_{T}$. The nondimensionalized Boussinesq equations are

$$\begin{eqnarray}&&\displaystyle \frac{1}{\Pr }\left(\displaystyle \frac{\partial \widehat{{\boldsymbol{u}}}}{\partial t}+\widehat{{\boldsymbol{u}}}\cdot {\rm{\nabla }}\widehat{{\boldsymbol{u}}}\right)=-{\rm{\nabla }}\widehat{p}+(\widehat{T}-\widehat{\mu }){{\boldsymbol{e}}}_{z}+{{\rm{\nabla }}}^{2}\widehat{u},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \widehat{T}}{\partial t}+\widehat{{\boldsymbol{u}}}\cdot {\rm{\nabla }}\widehat{T}+w({{se}}^{-t/{\tau }_{\mathrm{cool}}}+1)={{\rm{\nabla }}}^{2}\widehat{T},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \widehat{\mu }}{\partial t}+\widehat{{\boldsymbol{u}}}\cdot {\rm{\nabla }}\widehat{\mu }=\tau {{\rm{\nabla }}}^{2}\widehat{\mu },\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot \widehat{{\boldsymbol{u}}}=0,\end{eqnarray} \tag{ 16 }$$

where $s={T}_{0z}(t=0)/| {T}_{\mathrm{ad},z}|$. Note that the vertical potential density gradient is given by

$$\begin{eqnarray}&&\displaystyle \frac{d\widehat{\rho }}{{dz}}-\displaystyle \frac{d{\rho }_{\mathrm{ad}}}{{dz}}=-\left({{se}}^{-t/{\tau }_{\mathrm{cool}}}+1+\displaystyle \frac{d\widehat{T}}{{dz}}\right)+\displaystyle \frac{d\widehat{\mu }}{{dz}},\end{eqnarray} \tag{ 17 }$$

where $\tfrac{d{\rho }_{\mathrm{ad}}}{{dz}}=1$ in the system of units selected. The Ledoux threshold for convective instability then conveniently reduces to only requiring that the potential density gradient be zero.

In what follows, we solve these equations using the pseudo-spectral, double-diffusive convection code PADDI. This code is described in Traxler et al. (2011), and has been used extensively both in the geophysical (Stellmach et al. 2011; Garaud et al. 2015; Reali et al. 2017; Brown & Radko 2019) and astrophysical (Rosenblum et al. 2011; Mirouh et al. 2012; Brown et al. 2013; Garaud & Kulenthirarajah 2016, and many others) literature to date. PADDI uses periodic boundary conditions in all directions. As such, the initial conditions are constrained to be periodic as well, and must be chosen carefully to avoid any effect of the boundary conditions on the dynamics of the convective region that develops near the CO/ONe core–mantle interface. We have verified a posteriori for each run presented here, that the domain boundaries are far from the convective region and have no influence over it (see Figure 3, for instance).

Zoom In Zoom Out Reset image size

We select an initial mean molecular weight perturbation profile $\widehat{\mu }$ that is zero in the lower part of the domain (so the total mean molecular weight tends to ${\mu }_{{\rm{m}}}$, which is that of the CO core), and increases relatively sharply outward at some specified height, to model the presence of the CO/ONe core–mantle interface. Close to the upper boundary, the mean molecular weight perturbation is assumed to drop to zero again to satisfy the periodicity requirement. This upper transition does not affect the dynamics of the system much, as long as it is sufficiently far from the core–mantle boundary (see below for more on this topic). Since the background stable temperature gradient gradually decreases with time, our initial conditions for $\widehat{\mu }$ should also be selected to model a state that is not initially unstable to overturning convection, but that becomes so as the simulation proceeds.

To this end, our initial condition for the mean molecular weight is

$$\begin{eqnarray}&&{\widehat{\mu }}_{\mathrm{init}}(z)=\displaystyle \frac{{\rm{\Delta }}\mu }{2}\left[\tanh \left(\displaystyle \frac{z-{z}_{I}}{{\sigma }_{I}}\right)-\tanh \left(\displaystyle \frac{z-{z}_{T}}{{\sigma }_{T}}\right)\right].\end{eqnarray} \tag{ 18 }$$

The first $\tanh$ in ${\widehat{\mu }}_{\mathrm{init}}$ represents the core–mantle interface centered around z_I, with an initial destabilizing composition gradient

$$\begin{eqnarray}&&\displaystyle \frac{d{\widehat{\mu }}_{\mathrm{init}}}{{dz}}(z={z}_{I})=\displaystyle \frac{{\rm{\Delta }}\mu }{2{\sigma }_{I}}.\end{eqnarray} \tag{ 19 }$$

The size ${\rm{\Delta }}\mu$ and width σ_I of the mean molecular weight jump are two of the input parameters of the system, and must be carefully selected to achieve the desired parameter regime (see Section 3.3). The second $\tanh$ in ${\widehat{\mu }}_{\mathrm{init}}$ is required so that we satisfy the periodic boundary conditions required by our code. It is centered a little below the top of the computational domain, and the transition is selected to be very sharp to offer a strongly stabilizing buffer that prevents any remaining fluid motions in that region from moving through the top boundary. Our domain size will be selected to be sufficiently large, such that the turbulent motions of interest do not reach this transition region.

The initial conditions for the velocity fields are selected to be zero. A small amount of noise⁵ is added to the initial temperature and mean molecular weight profiles to drive the instability. We selected the initial amplitude of the random perturbations in the mean molecular weight to be 10⁻³, and in the temperature to be 10⁻¹. This is sufficiently small enough that the initial development of the instability is in the exponentially growing, linear phase. It is sufficiently large enough that we reach the more interesting nonlinear stage of the simulation on a timescale much shorter than that over which the gradients imposed in the initial condition would diffuse away. With these conditions met, the results will be independent of the initial choice of perturbation amplitude.

The nondimensional equations and initial conditions presented earlier contain 11 parameters, including some related to (1) the geometry of the system, characterized by the global domain dimensions $({L}_{x},{L}_{z})$ and the respective heights and widths of the transitions in the initial mean molecular weight profile ${z}_{I},{z}_{T},{\sigma }_{I}$ and σ_T, (2) the diffusion coefficients, characterized by the Prandtl number Pr and the diffusivity ratio τ, and (3) the system timescales, characterized directly by the global cooling time τ_cool, and indirectly by the initial temperature and mean molecular weight stratifications s and Δμ. These must therefore be carefully selected in order to achieve dynamics that resemble the ones expected to take place in the WD, while being mindful of computational restrictions on the simulation resolution and total integration time.

Starting with the diffusion coefficients Pr and τ, it is informative to determine what the latter are in the WD at the core–mantle interface. Characteristic temperatures and densities at this location, just prior to the onset of convection, are $T\approx 2\times {10}^{8}\,{\rm{K}}$ and $\rho \approx 1\times {10}^{7}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ (Brooks et al. 2017). We estimate the transport properties by evaluating existing literature estimates at these conditions for a 50/50 C/O mixture ($\bar{Z}\approx 7$, $\bar{A}\approx 14$). The electrons dominate the thermal conductivity and ${\kappa }_{T}\sim {10}^{3}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$ (Itoh et al. 1983). Likewise, the electronic contribution to the shear viscosity dominates and we have $\nu \sim {10}^{-1}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$ (Itoh et al. 1987). Thus, the Prandtl number of this material is $\Pr \sim {10}^{-4}$. From Beznogov & Yakovlev (2014), we get a diffusion coefficient ${\kappa }_{\mu }\sim 4\times {10}^{-3}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$. This implies the diffusivity ratio (reciprocal Lewis number) is $\tau \sim 3\times {10}^{-6}$.

We note that the initial size of the interface region is of the order of a pressure scale height $H\sim {10}^{8}\,\mathrm{cm}$. The magnitude of the destabilizing composition gradient can be expressed as the composition part of the Brunt–Väisälä frequency (i.e., Equation (1) assuming ${{\rm{\nabla }}}_{T}={{\rm{\nabla }}}_{\mathrm{ad}}$), which is ${N}^{2}\sim -0.1\,{{\rm{s}}}^{-2}$. This implies that the Rayleigh number is

$$\begin{eqnarray}&&\mathrm{Ra}=\displaystyle \frac{| {N}^{2}| {H}^{4}}{\kappa \nu }\sim {10}^{29},\end{eqnarray} \tag{ 20 }$$

where we used the diffusivities from the previous paragraph. Given this large Rayleigh number, convection is expected to be extremely efficient, which suggests that temperature should be fully mixed rather than remaining at the radiative gradient.

Unfortunately, such small values of Pr and τ are not computationally achievable in direct numerical simulations. Instead, we select for simplicity values that are closer to one (but still significantly below one). This will result in numerical simulations that qualitatively reproduce the correct dynamics, but cannot be expected to be quantitatively accurate.⁶ In order to quantify the impact of varying the parameters toward their true stellar values, we will compare the results of two sets of parameters: $\Pr =\tau =0.1$ and $\Pr =\tau =0.01$, respectively.

In the fingering regime, typical nondimensional turbulent structures have a size of order 10 (Garaud 2018), while in the convective regime we expect them to be as large as the convective zone. With the selection of Pr and τ made above, the size of the viscous and compositional boundary layers appearing at the edges of these turbulent structures will be of the order 0.1–1, so this sets the minimum nondimensional lengthscale that needs to be resolved. The total domain size will therefore be limited by the computational power available, given the required resolution. In what follows, we take L_x = 800 and L_z = 1200. We have verified that the horizontal size is large enough to contain many convective eddies (so horizontally averaged quantities are meaningful), and that the vertical size is large enough to have negligible influence on the system dynamics. With this vertical size, the convective region does not reach the domain boundaries or z_T during the simulation.

The selection of the remaining parameters is guided by the following considerations. We first need to ensure that the cooling time is substantially shorter than the fingering mixing timescale across the core–mantle interface. If this were not the case, it would contradict the idea that the thermohaline mixing is sufficiently inefficient, so as not to interfere with the flame (Denissenkov et al. 2013). In our nondimensional model, the turbulent diffusivity due to fingering convection is equal to $\tau {\mathrm{Nu}}_{\mu }$, where ${\mathrm{Nu}}_{\mu }$ is the Nusselt number associated with compositional transport in fingering convection, which has a typical value (Garaud 2018) of ${\mathrm{Nu}}_{\mu }\sim 10$ at the selected values of $\Pr$ and τ (except very close to the onset of overturning convection where it could be much larger), so we require

$$\begin{eqnarray}&&\displaystyle \frac{{\sigma }_{I}^{2}}{10\tau }\gg {\tau }_{\mathrm{cool}}.\end{eqnarray} \tag{ 21 }$$

We also have to ensure that the cooling time is much larger than the typical convective growth and turnover time, and this can be verified a posteriori. Note, however, that in the selected unit system, the Brunt–Väisälä frequency is proportional to ${\Pr }^{1/2}$, so the convective turnover timescale is expected to scale with ${\Pr }^{-1/2}$. As such, we must select a larger ${\tau }_{\mathrm{cool}}$ for the lower $\Pr$ simulations.

We also need to ensure that the initial condition is stable against overturning convection. The Ledoux criterion for stability, when expressed nondimensionally at $z={z}_{I}$ (which is the most unstable location in the domain at t = 0), implies

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}\mu }{2{\sigma }_{I}}\lt s+1.\end{eqnarray} \tag{ 22 }$$

Finally, the selected value of s needs to remain of the order of unity to realistically model the conditions within the WD, and the values of ${\sigma }_{I}$ and ${\sigma }_{T}$ must be much smaller than L_z.

Based on these considerations, we choose the parameters ${\rm{\Delta }}\mu =400$, ${\sigma }_{I}=40$, and s = 5. In order to ensure that the cooling timescale is much larger than the convective timescale, we pick ${\tau }_{\mathrm{cool}}=40$ when $\Pr =\tau =0.1$, and ${\tau }_{\mathrm{cool}}=120$ when $\Pr =\tau =0.01$. The remaining parameters are selected to be z_I = 500, z_T = 1100, and ${\sigma }_{T}=10$. The profile ${\widehat{\mu }}_{\mathrm{init}}(z)$, and its corresponding potential density gradient profile $d{\widehat{\rho }}_{\mathrm{init}}/{dz}\,-d{\rho }_{\mathrm{ad}}/{dz}=-(s+1)+d{\widehat{\mu }}_{\mathrm{init}}/{dz}$ at t = 0, are shown in Figure 3.

As required, the potential density gradient is negative everywhere at t = 0, so the system is stable against overturning convection (Ledoux-stable). Note that parts of the domain close to z_I are initially unstable to fingering convection, so we expect the fingering instability to first develop in these regions, followed by overturning convection once the stabilizing background temperature gradient has dropped below some critical threshold so that $d\bar{\rho }/{dz}-d{\rho }_{\mathrm{ad}}/{dz}$ becomes positive somewhere in the domain.

We integrated the set of Equations (13)–(16) with the initial conditions described in Section 3.2 for 12 cooling times, using a resolution of 1024 × 1536 Fourier modes (equivalently 3072 × 4608 grid points) for the $\Pr =\tau =0.1$ run, and for 6 cooling times using 2048 × 3072 Fourier modes (equivalently 6144 × 9216 grid points) for the $\Pr =\tau =0.01$ run. With the latter being a much more computationally demanding case, we were not able to run it for as long, relative to the cooling time. As we shall demonstrate, the qualitative evolution of the system is very similar in both cases when appropriately scaled, showing that the input parameters were correctly selected to be in the desired WD regime (e.g., where the cooling timescale is much faster than the fingering timescale, but much slower than the convective timescale). We have verified, by inspection of the Fourier spectrum (see the Appendix), that all scales of the simulations are fully resolved, from the global convective scale down to the viscous and compositional dissipation scales.

Figure 4 shows representative snapshots of the mean molecular weight perturbation field $\widehat{\mu }$ in the $\Pr =\tau =0.01$ run, when convection first sets in, and after one, three, and five cooling times. Profiles of the horizontally averaged potential density gradient $d\bar{\rho }/{dz}-d{\rho }_{\mathrm{ad}}/{dz}$, mean molecular weight $\bar{\mu }$, and potential temperature $\bar{T}-{T}_{\mathrm{ad}}=\bar{T}+z$ at these times are shown in Figure 3. The first snapshot shows the end of the fingering-only phase, and confirms that this instability is too weak to cause much vertical compositional mixing (as expected from the parameters selected) prior to the onset of overturning convection. It does imprint a particular initial size-scale on the developing convective eddies, however. That initial scale rapidly coarsens in the convective phase (second and third snapshots), consistent with the notion that the horizontal size of convective eddies is commensurate with the vertical extent of the convective zone (in the Boussinesq limit at least), which itself grows with time as the stabilizing background temperature gradient decays (see dashed black lines at one and three cooling times). Once fully developed, the turbulence exhibits a wide range of scales, as expected, from the double-diffusive scale ($l\sim 1$) to the size of the convective region. Inspection of the potential density gradient, potential temperature, and mean molecular weight profiles in Figure 3 shows that the convective zone is adiabatic (with $d\bar{\rho }/{dz}-1\simeq 0$), with substantial "noise" due to the fact that the simulations are only two-dimensional. The potential temperature and chemical composition of the convective zone is, however, not perfectly mixed, and small but finite compensating gradients of both quantities remain. This is interesting from a hydrodynamical point of view and deserves future investigations.

Zoom In Zoom Out Reset image size

The convective region stops growing after about four cooling times (with these selected parameters), at which point the background radiative temperature gradient ${{se}}^{-t/{\tau }_{\mathrm{cool}}}$ has become negligible, compared with the adiabatic temperature gradient (which is 1 in these units). The size of the convective region at that point is about 370 (in the selected units), and its boundaries are far from the domain boundaries, showing that the reason the convective region stops growing is not because it is computationally constrained to do so, but instead, because it has reached its equilibrium size. The growth of the size of the convective zone ${H}_{\mathrm{cz}}(t)$ is more easily visualized in Figure 5, together with that of the $\Pr =\tau =0.1$ simulation. To compute ${H}_{\mathrm{cz}}(t)$, we differentiate the horizontally averaged density profile $\bar{\rho }(z)$ with respect to z, compute the corresponding potential density gradient, then smooth the result using a 50 grid-point wide moving window. We then identify the locations of the points where the potential density gradient first crosses 0 from the bottom and from the top as the base and top of the convection zone, respectively.

Zoom In Zoom Out Reset image size

We see that the growth of ${H}_{\mathrm{cz}}(t/{\tau }_{\mathrm{cool}})$ is very similar in both runs, confirming that the growth of the convective zone is independent of the properties of convection itself (as long as the convective turnover timescale is short enough compared with ${\tau }_{\mathrm{cool}}$), but instead, only depends on the global structure of the star. In fact, a good estimate for ${H}_{\mathrm{cz}}(t)$ can be obtained by requiring that the total jump in potential density across the convection zone be zero. Since the total jump in potential temperature is ${\rm{\Delta }}{T}_{\mathrm{cz}}(t)={H}_{\mathrm{cz}}(t)\left[s\exp (-t/{\tau }_{\mathrm{cool}})+1\right]$, and the total jump in mean molecular weight is ${\rm{\Delta }}{\mu }_{\mathrm{cz}}={\widehat{\mu }}_{\mathrm{init}}({z}_{I}+{H}_{\mathrm{cz}}/2)\,-{\widehat{\mu }}_{\mathrm{init}}({z}_{I}-{H}_{\mathrm{cz}}/2)$ (assuming that regions outside of the convective zone have not been mixed yet), then ${H}_{\mathrm{cz}}(t)$ can be found by solving the equation

$$\begin{eqnarray}&&{H}_{\mathrm{cz}}(t)\left[s\exp (-t/{\tau }_{\mathrm{cool}})+1\right]\simeq {\rm{\Delta }}\mu \tanh \left(\displaystyle \frac{{H}_{\mathrm{cz}}(t)}{2{\sigma }_{I}}\right),\end{eqnarray} \tag{ 23 }$$

using a Newton–Raphson relaxation algorithm. The solution is also shown in Figure 5. We see that this provides a relatively good approximation to the extent of the convective zone measured in both simulations. Deviations from the theoretical predictions can be attributed to (1) uncertainties in the measurement of ${H}_{\mathrm{cz}}$ due to the large fluctuations in $d\bar{\rho }/{dz}$, and (2) to mixing beyond the edge of the convection zone due to overshoot or fingering, which is ignored in this model.

Figure 6 shows another more quantitative comparison of the $\Pr =\tau =0.1$ and $\Pr =\tau =0.01$ simulations, by looking at the evolution of the total volume-averaged rms velocity ${u}_{\mathrm{rms}}$ as a function of time. We see that the two are very close to one another when time is scaled by the global cooling time, and when ${u}_{\mathrm{rms}}$ is scaled by the Prandtl number, at least up to $t/{\tau }_{\mathrm{cool}}\simeq 6$. We also see that ${u}_{\mathrm{rms}}$ grows almost linearly with time at early times, mirroring the growth of ${H}_{\mathrm{cz}}(t)$. This confirms that the typical convective velocity in the system is proportional to ${H}_{\mathrm{cz}}N$ where, as discussed earlier, the Brunt–Väisälä frequency $N\propto {\Pr }^{1/2}$.

Zoom In Zoom Out Reset image size

As the convection zone slowly equilibrates, we see in Figure 4 the development of extended regions on either side that are mixed by fingering convection, consistent with the predictions of Brooks et al. (2017). As these regions expand in size, they mix the areas immediately adjacent to the convection zone and the latter shrinks in size (see Figure 5, in particular the $\Pr =\tau =0.1$ case). This causes the volume-averaged rms velocity in the convection zone to decrease accordingly (see Figure 6). Since turbulent mixing by fingering convection at $\Pr =\tau =0.1$ is more efficient than at $\Pr =\tau =0.01$, the temperature and compositional gradients that form on either side of the convection zone are shallower in the former case, and steeper in the latter (see Figure 3). As a result, the decrease in ${H}_{\mathrm{cz}}$ and ${u}_{\mathrm{rms}}$ is faster in the $\Pr =\tau =0.1$ case than in the $\Pr =\tau =0.01$ case. In other words, once cooling is no longer actively driving the convection, the time-evolution of the system becomes controlled by fingering-induced mixing at the edge of the convection zone instead. A similar transition also occurs in the MESA models of Brooks et al. (2017), as the initially dominant mixing from overturning convection gradually fades and thermohaline mixing takes over.