Fast and Slow Crystallization-driven Convection in White Dwarfs

We investigate crystallization-driven convection in carbon–oxygen white dwarfs. We present a version of the mixing length theory that self-consistently includes the effects of thermal diffusion and composition gradients, and provides solutions for the convective parameters based on the local heat and composition fluxes. Our formulation smoothly transitions between the regimes of fast adiabatic convection at large Peclet number and slow thermohaline convection at low Peclet number. It also allows for both thermally driven and compositionally driven convection, including correctly accounting for the direction of heat transport for compositionally driven convection in a thermally stable background. We use the MESA stellar evolution code to calculate the composition and heat fluxes during crystallization in different models of cooling white dwarfs, and determine the regime of convection and the convective velocity. We find that convection occurs in the regime of slow thermohaline convection during most of the cooling history of the star. However, at the onset of crystallization, the composition flux is large enough to drive fast overturning convection for a short time (∼10 Myr). We estimate the convective velocities in both of these phases and discuss the implications for explaining observed white dwarf magnetic fields with crystallization-driven dynamos.


Introduction
As white dwarfs (WDs) cool over time, a first-order phase transition leads to the solidification of their degenerate interiors (van Horn 1968).This crystallization process occurs from the core outward.In carbon-oxygen WDs, the transition from liquid to solid preferentially retains oxygen in the solid core, creating an unstable, buoyant region on top of it as lower-density carbon is released.Thus, convection is induced as the crystallization front grows outward, changing the composition profile (Stevenson 1980;Mochkovitch 1983;Isern et al. 1997;Fuentes et al. 2023).The latent heat released during the phase transition can significantly slow down the cooling of the star, and the delay is increased even more by the gravitational energy released as the elements are redistributed.The cooling delay leads to an overdensity in the Hertzsprung-Russell diagram, known as the Q-branch, that has been revealed by Gaia (Tremblay et al. 2019), although an additional source of cooling delay seems to be required to explain the observed overdensity (Cheng et al. 2019).
The crystallization process has also been linked to magnetism in white dwarfs.It has been suggested that compositionally driven convection during crystallization can drive a dynamo (e.g., Isern et al. 2017;Ginzburg et al. 2022).This would explain why most magnetic white dwarfs appear at low temperatures and luminosities (e.g., Liebert et al. 2003;Sion et al. 2014) as well as their incidence in different close binary systems (Belloni et al. 2021;Schreiber et al. 2021Schreiber et al. , 2022)).However, whether convection is efficient enough to explain the large intensity of the observed magnetic fields ( obs ∼ 0.1-1000 MG) is still under debate.If convection is efficient and the gravitational energy released during crystallization is used to power the magnetic field, Isern et al. (2017) showed that the saturated dynamo scalings of Christensen et al. (2009) (see also Christensen 2010) can explain field strengths in the MG range.Taking rotation into account could potentially explain larger fields (Isern et al. 2017;Ginzburg et al. 2022).
Although the energetics of a crystallization-driven dynamo seem reasonable, it is not clear whether the convective velocity is large enough.For a saturated dynamo, we expect  2  ∼  2 /4, which gives   ∼ 100 cm s −1 for  ∼ 1 MG (Ginzburg et al. 2022).However, this velocity is much larger than expected based on estimates from mixing length theory (MLT) applied to white dwarfs, in which the large thermal conductivity reduces the convective efficiency by reducing the effective buoyancy of rising fluid elements (Mochkovitch 1983;Isern et al. 1997).The kinetic energy flux is then a small fraction of the total energy flux.Mochkovitch (1983) estimated convective velocities in the range 10 −6 -10 −1 cm s −1 , depending on the star's rotation rate.For such small velocities, if the kinetic energy carried by convection powers the magnetic field, the intensity of the field is restricted to a few Gauss, much lower than observed.
Using both MLT and hydrodynamic simulations, we recently showed that there are two modes of compositionally driven convection, depending on the magnitude of the composition flux (Fuentes et al. 2023, hereafter F23).When the composition flux is larger than some threshold, convection occurs in a fast mode in which the convective motions are too rapid for thermal diffusion to be important.For this fast mode of convection, the velocity estimates are similar to those of Ginzburg et al. (2022).Below the threshold composition flux, convection is slow, with thermal diffusion acting to reduce the convective velocities to the small values found by Mochkovitch (1983).As a first estimate, F23 used cooling models of white dwarfs to obtain the composition flux, finding that, except for a brief period at the onset of crystallization, convection in white dwarfs occurs in the slow mode, with correspondingly small convective velocities.However, these estimates were approximate since the models used did not include phase separation self-consistently, and F23 assumed that the inward convective and outward conductive heat fluxes exactly balance (zero net heat flux).
Current evolution models of cooling WDs that include crystallization do not explicitly follow convection.Instead, it is usually assumed that the liquid phase above the core is fully mixed, because the relevant timescales considered for the evolution are much larger than the convective timescale (see, e.g., Althaus et al. 2012for LPCODE, Salaris et al. 2022 for BaSTI, Bauer 2023 for MESA).The energy contribution from the redistribution of chemical elements is then added to the structure equations as a source term following the prescriptions in Isern et al. (1997Isern et al. ( , 2000)).The energy released by the phase transition itself, i.e., the latent heat, is naturally included as a source term in the local solidification front.Such a treatment simplifies the computation of the cooling while still including the effect of phase transition on the total energy released and therefore cooling delay.However, the assumption that the layers above the solid core are instantaneously mixed means that the parameters of convection are not calculated.
One might imagine that, to follow convection in a cooling white dwarf, we could drop the assumption of instantaneous mixing and instead allow the convective mixing routine in the code to do the mixing for us.However, the standard prescriptions for convection in stellar evolution codes are not applicable to compositionally driven convection in white dwarfs.A prescription for inefficient convection is usually included in stellar evolution codes in which MLT is modified to include radiative losses (e.g., see Kippenhahn et al. 2013 which is based on Böhm-Vitense 1958).However, this applies only for overturning convection in which the temperature gradient is superadiabatic (for example, in surface convection zones of stars).In the slow regime of convection in white dwarfs, thermal diffusion is efficient enough that the temperature gradient remains subadiabatic.This is the regime of thermohaline convection (e.g., Kippenhahn et al. 1980Kippenhahn et al. , 2013)), which is typically modeled with a separate prescription that transports composition but not heat (e.g., Kato 1966;Kippenhahn et al. 1980;Langer et al. 1983).Even in the fast regime of compositionally driven convection where thermal losses can be neglected, standard MLT implementations do not apply, since they usually assume that the convective heat transport is outward, whereas compositionally driven convection in white dwarfs transports heat inward (F23).
In this paper, we discuss the formalism needed to follow convection and the associated composition and heat transport in 1D cooling white dwarf models.We extend the MLT developed by F23, which assumed zero net heat flux, to arbitrary heat and composition fluxes.We show that the cubic equation that is usually solved to obtain the convective velocity in MLT is replaced by a fifth-order equation when composition gradients and thermal losses are correctly accounted for.A similar equation was derived by Umezu & Nakakita (1988) for thermally driven overturning convection in massive stars; here, we extend this result to also include compositionally driven overturning convection and thermohaline convection.
We start in Section 2 with a discussion of the standard implementations of MLT for inefficient convection and its implementation in stellar evolution codes and then present our formulation for studying compositionally driven convection.In Section 3 we explore the regimes of compositionally driven convection as a function of the total heat flux and composition flux.In Section 4, we discuss the application of our theory to cooling models of WDs.In particular, we find that convection occurs in the fast regime for a short time after the onset of crystallization (≲ 10 Myr), but then quickly transitions and remains in the slow regime for the rest of the evolution.In Section 5, we conclude and discuss the implications for future white dwarf models and crystallization-driven dynamos.

Mixing Length Theory
In this section, we first review the standard treatment of inefficient convection in the framework of MLT (Section 2.1), and then discuss the ways in which it has been extended to include composition gradients in different stellar evolution codes (Section 2.2).Finally, we develop a general MLT that covers the regimes of thermohaline and overturning convection that we can apply to white dwarf interiors (Section 2.3).Additionally, we show how our MLT prescription relates to previous implementations of thermohaline convection (Section 2.4)

Inefficient Convection without Composition Gradients
We first discuss the usual implementation of inefficient convection without composition gradients, following the prescription of Cox & Giuli (1968).
The density contrast of convecting fluid elements is written as where we have used the notation of Kippenhahn et al. (2013) where  is the local gravitational acceleration.The effect of thermal diffusion is to set ∇  .If heat loss is negligible, fluid elements move adiabatically and ∇  = ∇ ad ; if heat loss is very efficient, the fluid element experiences the same gradient as the background, ∇  → ∇.The value of ∇  between these limits is determined by how quickly the fluid element moves.Considering the exchange of energy between the fluid element and the surroundings gives (Kippenhahn et al. 2013 based on Böhm-Vitense 1958) where  0 is a geometric term that depends on the assumptions for the shape of the fluid element, and is the thermal diffusivity, where  is the opacity,  the radiation constant, and  the speed of light.The dimensionless quantity Γ is usually taken as a measurement of the convective efficiency since it describes how closely convecting fluid elements follow the background gradient and therefore how easily they are accelerated by the buoyancy forces (Jermyn et al. 2022).Note that Γ is proportional to the Péclet number Pe ≡   ℓ/  (see Eq. (4) of F23) which measures the ratio of thermal diffusion timescale to convective turnover time.The adiabatic limit (efficient convection) corresponds to Γ ≫ 1, while inefficient convection occurs when Γ ≪ 1.Using Equation (2) for the velocity gives where we define (5) Note that, for a constant value of , we have The convective heat flux is written as   =    =    , where we relate the heat and entropy excess using the second law of thermodynamics  =  .Also, as long as the fluid is compositionally homogeneous,  =   /, where   is the specific heat capacity at constant pressure, giving The radiative heat flux is given in terms of the temperature gradient as It is also useful to define the gradient ∇ rad that measures the total heat flux, Using Equation (2) for the velocity, we obtain a relation between the temperature gradients Thus, when there are no composition gradients, the convective heat flux is proportional to (∇ − ∇  ) 3/2 .If the total heat flux (given by ∇ rad ) and the adiabatic gradient ∇ ad are known quantities, then ∇, ∇  , and Γ can be found using Equations (3), ( 4) and (9) (given the constants  0 and ).The solution for Γ is the only real and positive root of the following third-degree polynomial equation2 Denoting the root as Γ 0 , the solutions for the temperature gradients are where is another measure of the convective efficiency.
MLT is typically implemented in stellar evolution codes by solving Equation (10) or its equivalent to determine the temperature gradient ∇ given the total heat flux ∇ rad (e.g., the mlt routine in MESA, Paxton et al. 2011).When  → 1, convection is efficient with Γ ≫ 1 and the temperature gradient is close to adiabatic ∇ → ∇ ad .When  → 0, convection is inefficient, Γ ≪ 1, and radiation carries most of the heat with a superadiabatic gradient ∇ → ∇ rad .

Implementation of Composition Gradients in Stellar
Evolution Codes Composition gradients modify the criterion for the onset of convection from the Schwarzschild criterion ∇ > ∇ ad to the Ledoux criterion where =1    ∇   represents the composition gradients3,    =  ln / ln   | , , ≠ , and ∇   =  ln   / ln | ★ for  species with mass fractions   .The mass abundances must fulfill  =1   = 1, which allows us to eliminate one of the abundances in the ∇ com definition, whose summation has been written from 1 to −1.Depending on the sign of ∇ com , composition gradients can be either stabilizing (∇ com < 0) or destabilizing (∇ com > 0).
When the temperature gradient needed to carry the heat flux by radiation ∇ rad is unstable according to the Ledoux criterion, convection will occur and MLT can be applied.Different stellar evolution codes have taken different approaches to incorporating the term ∇ com in the basic Equations ( 3), (4), and (9).For instance, in MESA, the definition of Γ in Equation ( 3) is modified by replacing the adiabatic gradient with ∇ L , while the expressions for convective velocity and flux (Eqs.( 4) and ( 9)) are left unchanged (see Jermyn et al. 2023 Equations ( 18), ( 19) and ( 25) for flux, Γ, and   respectively).This gives a convective flux ∝ (Γ/(Γ + 1)) 3/2 (∇ − ∇  ) 3/2 .
As well as being different from each other, these implementations have a number of problems when applied to compositionally driven convection in white dwarfs.First, as we discuss in Appendix A, the composition terms in the entropy excess, which determines the heat flux, are a small correction, so that writing the heat flux in terms of the enthalpy excess overestimates the effect of the composition terms.Second, these implementations do not correctly handle the case where convection is driven by large destabilizing composition gradients that overcome a stable thermal gradient.When   ∝ (∇ − ∇  + ∇ com ) 3/2 , the convective heat flux is always outward when the Ledoux criterion is satisfied, whereas convection should transport heat inward in a thermally stable background (F23).This implies that the equations being used are not internally self-consistent; indeed, we will see below that the cubic equation that is usually solved to obtain Γ (e.g.eq.(A15) of Vazan et al. 2015) should be replaced by a quartic or quintic equation when composition gradients are included.
When the Ledoux criterion is not satisfied (∇ < ∇  ), the fluid is overall stable to convection and MLT gives Γ = 0 (no convection).However, double-diffusive instabilities can still drive mixing if either the thermal gradient or the com-3 In this notation, ∇ com is equivalent to − in Equation ( 6) of Paxton et al. (2013), where the implementation of composition gradients in MESA is described.This quantity is also usually written in terms of the mean molecular weight , which simplifies the summation term, ∇ com = − ( /  ) ∇  (see Maeder 2009;Kippenhahn et al. 2013), where  =  ln / ln  | , and  = − ln / ln  | ,  .position gradient is unstable.If the composition gradient is overall stabilizing, but the thermal gradient is unstable, this leads to semiconvection.In the case relevant for white dwarf crystallization, we have a stabilizing thermal gradient, but an unstable composition gradient, leading to thermohaline convection.These different kinds of double-diffusive convection are usually handled with a mixing prescription that is separate from MLT; see, e.g., Langer et al. (1983) for semiconvection or Kippenhahn et al. (1980), Brown et al. (2013), andLattanzio et al. (2015) for thermohaline convection.In particular, for thermohaline convection, the mixing is usually done by defining an effective diffusion coefficient for composition and it is assumed that all the heat transport is by radiative diffusion.Therefore, in regions undergoing thermohaline convection, there is no consideration of the convective heat transport.In the next section, we present an extension of MLT that smoothly transitions between the thermohaline and overturning convection regimes.

MLT Including Compositionally Driven Convection
Having motivated the need for a more complete version of mixing length theory to model compositionally driven convection, we now extend the MLT derived by F23 to arbitrary heat and composition fluxes.Whereas F23 assumed that the inward convective heat flux and outward conductive flux were always in balance, our goal here is to generalize this so that, given the heat flux and composition flux, we can predict the temperature and composition gradients required for the stellar model.
We proceed as in Section 2.1, but including composition gradients in each step.The density contrast in the convection zone is Since the buoyant acceleration is proportional to − /, continued convection requires   < 0 or (16) In the limit of rapid convection where thermal diffusion can be neglected, ∇  → ∇ ad , and condition ( 16) is just the usual form of the Ledoux criterion (Equation ( 14)).In the case of slow thermohaline convection, however, thermal diffusion reduces ∇  to a value well below ∇ ad , so that an unstable composition gradient ∇ com > 0 can sustain convection even though the adiabatic Ledoux criterion is not satisfied.
The convective velocity is differing from Equation (2) due to the addition of ∇ com .Using the definition of Γ from Equation (3), we obtain the general form of the convective efficiency which replaces Equation (4).Equation (3) still applies even with composition gradients, since it can be obtained using just the temperature excesses (Maeder 2009;Kippenhahn et al. 2013).As mentioned before, the contribution of the composition term to the entropy excess can be neglected in the case of WD crystallization (Appendix A), so that Equation ( 6) for the convective heat flux applies in this case as well.Using Equation ( 17) for the convective velocity, we find which replaces Equation ( 9).Note that from Equation ( 19), compositionally driven convection in a subadiabatic thermal background (∇ < ∇  ) has a heat transport directed inward (as in F23).
Combining Equations ( 3), ( 18), and ( 19), we find the fourthdegree equation which reduces to Equation ( 10) when ∇ com = 0. Unlike Equation ( 10), Equation ( 20) does not always have a unique, real, and positive root for given values of ∇ rad , ∇ ad , and ∇ com .However, even though composition gradients are usually incorporated into MLT by using ∇ com , a more natural formulation is to use the composition flux   =   ∇  (ℓ/2  ) as an input parameter rather than composition gradient4.The reason for this is that, in Henyey-type codes, the heat and composition fluxes are fundamental variables in the energy and mass conservation equations, so we generally need to calculate the temperature and composition gradients given the fluxes and not the other way around.In the case of white dwarf crystallization, the light element flux is set by the cooling rate of the core.F23 showed that a natural dimensionless parameter to use that measures the composition flux is The composition flux that gives  = 1 has an associated convective heat flux that is equal to the conductive heat flux down the adiabat.Thus,  gives a measure of how much the convection tries to adjust the temperature gradient.For their case of zero net heat flux (∇ rad = 0), F23 found a unique, positive and real solution for Γ as a function of , with a rapid transition between the regimes of slow convection with Γ ≪ 1 and fast convection with Γ ≫ 1 occurring at  ∼ 1.
Rewriting Equation (20) in terms of  rather than ∇ com gives 4 Here, for WD crystallization, we consider a two-species mixture, with abundances  of carbon and  = 1 −  of oxygen; the composition gradient is then which generalizes the results of F23 to give Γ as a function of the two fluxes ∇ rad and  (setting ∇ rad = 0, Equation ( 22) reduces to Equation (25) of F23).Similar quartic and quintic equations were derived by Umezu & Nakakita (1988) for overturning convection in massive stars (compare their Equations ( 25) and ( 28)).It is straightforward to verify that the same relations between the gradients and Γ still hold as we did previously, so that once a solution for Γ is obtained, we can find the gradients ∇, ∇  , and ∇ com using Equations ( 11), (12), and ( 21).We explore the solutions to this equation in the next section.

Transition from Efficient to Inefficient Convection
In the following, we show that thermohaline convection can be understood as inefficient compositional convection (Γ ≪ 1).Recall that continued convection (either compositional overturning convection or thermohaline convection) demands condition ( 16) is satisfied.For the particular case of thermohaline convection, where ∇ − ∇ ad + ∇ com < 0 (Ledoux stable), Equation ( 16) requires that both the thermal and composition gradients must simultaneously fulfill Using Equation (3), the condition above can be rewritten in terms of the efficiency Since ∇ com > 0 and ∇ − ∇ ad < 0, Equation ( 24) requires the convective efficiency to be small Γ ≪ 1 (which is equivalent to having a thermal adjustment time  therm much smaller than the convective turnover time  conv ).This is consistent with the fact that the thermohaline instability develops in a thermally stable background, where thermal diffusion is important.Now, we show that our formalism gives a mixing coefficient for thermohaline convection that is of the same form as the prescriptions used for thermohaline convection in stellar evolution codes.In the limit Γ ≪ 1, the criterion for continued convection (Equation ( 16)) is marginally satisfied.Then, Equation ( 18) gives ∇ com ≈ (∇  − ∇), and from the middle term of Equation (3) we can write If we define the mixing coefficient as  =   ℓ and use the definition of Γ in the last term of Equation (3) we find where the ratio  0 = (∇ ad − ∇)/∇ com (as in Equation 43of Brown et al. 2013).This expression for  is the same as that -10 3 -10 0 -10 3 0 10 3 10 0 10 3 -10 2 -10 0 -10 2 -10 4 0 10 4 10 2 10 0 10 2 rad / ad
We make clear that our set of equations does not directly provide a solution for the mixing length ℓ.In Section 4, we use the numerical solutions of our MLT equations for Γ and the relevant gradients to explore the effect of different assumptions of ℓ in the inefficient regime.

The Different Regimes of Compositionally Driven Convection
We now use the formalism established in the previous section to investigate the different regimes of compositionally driven convection.We assume that the heat flux and composition flux are specified (as dimensionless parameters ∇ rad /∇ ad and  respectively; see Equations ( 8) and ( 21)) and calculate the convective efficiency Γ (related to the convective velocity and Peclet number as in Equation ( 3)) and the temperature and composition gradients ∇ and ∇ com .Figure 1 displays the numerical solutions of Equation ( 22) for the choice of parameters5  0 = 9/4,  = 2.5 × 10 4 , and ∇ ad = 1/3.Except for situations where the condition   < 0 (Equation ( 16)) was not fulfilled (i.e., regions where convection does not occur and Γ is set to 0), there is a unique solution for Γ(∇ rad , ) in almost the entire parameter space.There is a second solution in the top left region of Figure 1 ( < 0 and ∇ rad /∇ ad > 0), but it shows an unphysical discontinuity across  = 0 (see Appendix B), and so we do not show it here.
The different regimes of convection are labeled in Figure 1.Whereas we usually think of the boundaries between different convective regimes in terms of the gradients (for example, whether or not the gradients satisfy the Ledoux criterion), here we have parameterized the convection by the heat and composition fluxes.In this case, it is helpful to write Equation ( 16) for   < 0 in terms of the fluxes.Replacing ∇ − ∇  and ∇ com from Equations ( 12) and ( 21), and multiplying by  0 Γ we can write Equation ( 16) as follows: Finally, dividing by ∇ ad and using the definition of  from Eq. ( 13) to write Since  ranges from 0 to 1, we see that a total heat flux ∇ rad > ∇ ad or a positive composition flux  > 0 destabilize, whereas a total heat flux ∇ rad < ∇ ad or negative composition flux  < 0 stabilize.When either ∇ rad or  are large and positive, we have overturning convection (Ledoux unstable) with Γ ≫ 1.The red curve in Figure 1 bounds the region of thermohaline convection with Γ ≪ 1, which is unstable according to Equation ( 16) (  < 0) but Ledoux stable (∇ − ∇ ad + ∇ com < 0).To distinguish between each convective regime, we use Equations ( 11), (12), and (21) to obtain numerical solutions for ∇, ∇  , and ∇ com .For example, the red curve is defined by the solutions that satisfy ∇ − ∇ ad + ∇ com = 0, where the bounded zone contains all the solutions that fulfill the conditions for thermohaline convection, i.e., ∇ − ∇  + ∇ com > 0 > ∇ − ∇ ad + ∇ com .
At the top left, we see that overturning convection with ∇ rad > ∇ ad shuts off if there is a large enough stabilizing composition flux (negative ).This is the region in which semiconvection would occur (bordered by the magenta curve in Figure 1, which we draw similarly to the red curve for the thermohaline zone).Note, however, that as semiconvection is not included in our theory, Equation ( 16) gives   > 0 in this region and Γ = 0.When both  < 0 and ∇ rad < ∇ ad (lower left), convection shuts off completely and Γ = 0.
-10 3 -10 0 -10 3 0 10 3 10 0 10 3 -10 2 -10 0 -10 2 -10 4  Figure 1 shows that compositionally driven convection has two distinct regimes: slow thermohaline convection with Γ ≪ 1 and fast efficient overturning convection with Γ ≫ 1.As previously found by F23, the transition from slow to fast is extremely rapid and occurs at  ≈ 1 as long as the total heat flux is not too negative (∇ rad > −∇ ad ).If ∇ rad < −∇ ad , the total inward heat flux, which is stabilizing, can be large enough to compensate for the destabilizing outward composition flux, reducing Γ.This moves the slow-fast transition to larger values of  ≈ −∇ rad /∇ ad (lower right of Figure 1).
Figure 2 shows the behaviour of the gradients ∇ (left panel) and ∇ com (right panel) in the different regimes.In the stable regions where Γ = 0 (see Figure 1), ∇ is just equal to ∇ rad , meaning that the energy transport occurs only by radiation.In regions where Γ ≫ 1, ∇/∇ ad → 1, independently of the nature of the dominant destabilizing factor, i.e., outward composition or heat flux.For the thermohaline regime, the value of ∇ depends on both fluxes, increasing from ∇ rad to ∼ ∇ ad in the range of  ∼ 0-1.The transition becomes more abrupt as the stabilizing heat flux is more negative.In the thermohaline regime with large inward heat flux, ∇ changes sign.The transition from positive to negative can be understood by setting ∇ = 0 in Equation ( 11) and using the Γ ≪ 1 solution of Equation ( 22),  ≈  (1 − ∇ rad /∇ ad ), which gives  = −∇ rad /∇ ad for ∇ = 0.This agrees with the location of the change in sign of ∇ in the left panel of Figure 2.
The right panel of Figure 2 shows the solutions for ∇ com .Its sign is given by the direction of the composition flux , and its magnitude also is dominated by how large or small is .While ∇ com is generally an increasing function of , there is a sudden drop in its magnitude when convection transitions from slow to fast convection.Overall, the composition gradient depends mostly on  and is more or less independent of ∇ rad .The exception is for strong inward heat flux ∇ rad < −∇ ad , where the composition gradient increases as the inward heat flux becomes stronger at fixed .
Figure 3 shows the convective heat flux   , normalized by  ad =      ∇ ad /  , the conductive heat flux along the adiabat (this is equivalent to normalizing Equation ( 19) for ∇ rad − ∇ by ∇ ad ).For thermal overturning convection when ∇ rad > ∇ ad , the convective heat flux is outward, whereas for compositionally driven convection, the convective heat flux is inward.This follows from the fact that   ∝ (∇ − ∇  ), so it is positive when the regime is superadiabatic and negative for subadiabatic.For overturning convection,   depends mainly on ∇ rad , and for thermohaline convection,   depends mainly on .

Application to White Dwarf Crystallization
In F23, we used models of cooling white dwarfs from the MESA code to estimate the parameters of convection during the crystallization phase.However, those models did not include phase separation and the convection theory assumed that the conductive and convective heat fluxes exactly balanced, giving ∇ rad = 0.In this section, we improve on those estimates by using cooling models that include phase separation, following Bauer (2023), and using the total heat flux and composition flux extracted from the model to determine the convective regime and obtain the parameters of convection.We note that, in light of the mixing length theory presented in Section 2.3, a self-consistent calculation of  and Γ requires implementing such formalism into 1D evolution models.This will be accomplished in future work.
We used MESA version r23.05.1 (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019;;Jermyn et al. 2023) to run cooling mod-els6 of white dwarfs incorporating phase separation of the C/O core following Bauer (2023).We adopt two different composition profiles: (1) a realistic profile from the stellar evolution of a 0.6 ⊙ WD obtained directly from the settled model of the test suite wd_cool_0.6M,and (2) an initial profile from Bauer (2023) corresponding to a WD with a uniform core abundance of 50% carbon (C) and 50% oxygen (O).To simplify the calculation, we use a nuclear network containing only C, O, He, and H (with nuclear reactions turned off), so that the internal structure of the star is composed of a C/O core surrounded by a He/H atmosphere.The phase separation calculation takes into account the gravitational energy released due to the redistribution of elements (Bauer 2023) and the latent heat released during the phase transition is obtained from the Skye EOS in Jermyn et al. (2021).
To calculate the flux of light elements we first identified the location of the crystallization front in each model profile.Following the calculation of Bauer (2023), we use the phase value  = 0.9 from the Skye EOS (Jermyn et al. 2021) to set the solid core boundary ( = 0 represents the liquid phase, while  = 1 is the solid phase).
Then, following F23, we compute the composition flux as where   is the rate of growth of the core mass,   is the radius of the core, and Δ melt is the increase of C in the liquid phase relative to the solid, which is obtained from the liquidsolid composition difference given by the C/O phase diagram of Blouin et al. (2020) and Blouin & Daligault (2021).Once we estimate   , we compute the normalized composition flux  from the middle term in Equation ( 21), using the output profiles of the models to estimate the thermodynamic properties (, ,   ,   , ∇ ad ).To ensure that these parameters correspond to the liquid just outside the crystallization front in the convection zone, they were taken the model profiles at the location of a phase value  = 0.01.The local value of ∇ rad at  = 0.01 representing the total heat flux is also taken from the output profiles.We checked that varying this value of  in the range 10 −4 to 0.5 does not significantly change the results.
The results are shown in Figures 4 and 5 for two different white dwarf masses, 0.6 and 0.9  ⊙ , with the stellar evolution abundance profile.We also show a 0.6  ⊙ model with the 50/50 C/O profile to show the effect of changing the abundance profile.Figure 4 shows the evolution of the heat and composition fluxes in the -∇ rad /∇ ad plane.The color map shows the corresponding convective velocity   obtained from the efficiency parameter Γ using Equation (3), assuming  0 = 9/4 and taking typical values for , ∇ ad , and   along each trajectory.Figure 5 shows the time evolution of the solid core mass, composition at the liquid-solid interface, and the composition and heat fluxes.We find in all models, that for most of the evolution of the star, convection is in the slow thermohaline regime, where   ≲ 10 −5 cm s −1 , in agreement with F23.However, there is a short phase (≲ 10 Myr) of fast overturning convection ( > 10) at the beginning of crystal- lization, with convective velocities   ∼ 10 − 100 cm s −1 (see the inset in panel (c) of Fig. 5).The large composition flux at early times is due to the larger rate of crystallization in the first stages of the cooling (panel (a) of Figure 5) and the small core radius at that time.
The value of  at the early times following crystallization depends on the temporal resolution used.Smaller core values could be achieved if the initial crystallization time is solved more precisely.In this case, we used a time step of ∼ 10 5 yrs in the early stages for all the models, allowing us to solve an efficient regime that lasts only a few Myr.While we do not attempt to show how small the core could be, in a companion paper, Fuentes et al. (2024), we checked that using an analytical estimation for  also predicts an efficient convection regime that occurs during ∼ 10 Myr from the onset of crystallization.
Panel (d) of Figure 5 shows the total heat flux ∇ rad in the liquid just outside the crystallization front.The net heat flux is set by the combination of cooling luminosity, latent heat, and the gravitational energy released during the phase separation.At first, these contributions produce a sudden increase in the total outward heat flux at the onset of crystallization.In the 0.6  ⊙ white dwarf, the latent heat and gravitational energy from phase separation maintain the total heat flux at an approximately constant level as crystallization continues.As expected (Bauer 2023), the more massive 0.9  ⊙ white dwarf cools more quickly, and ∇ rad continues to decrease, although at a slower rate than if latent heat and phase separation were not included.
While ∇ rad /∇ ad moves in a tight range of values, keeping approximately constant, the composition flux  steadily decreases.At  ∼ 8 Gyr, we observe an interesting dependence on the abundance profile.For the stellar evolution C/O profile,  drops rapidly and becomes negative, whereas  becomes small but remains positive for the 50/50 C/O profile.This difference is a direct consequence of the C/O phase diagram.
According to the analytical fit presented in Blouin & Daligault (2021), when the number abundance of oxygen in the liquid falls below  ℓ O ≲ 0.18, the contrast between the oxygen abundance in the liquid and solid becomes negative (Δ O < 0).This means that the liquid phase is O-enriched and the solid is enhanced in C, creating a compositionally stable region at the base of the C/O mixture.This transition can be seen in panel (b) of Figure 5, where we indicate the normalized O mass fraction X O /(X C + X O ) ≈ 0.219 corresponding to the azeotrope in the phase diagram at  ℓ O ≈ 0.18.Comparing with panel (c), we see that  becomes negative when the O abundance drops below the azeotropic value.In contrast, the normalized mass fraction of oxygen in the models with 50/50 C/O initial abundance never crosses the limit X O /(X C + X O ) ≈ 0.219, and therefore  remains positive.This difference is related to the different composition profiles in the two cases.The stellar evolution profile has a central region with uniform abundances surrounded by an outer region in which the oxygen abundance drops (see Figure 8 of Bauer 2023).As the convection zone expands into this outer region, the oxygen abundance drops below the azeotropic value.In contrast, in the 50/50 C/O model, the flatter oxygen abundance profile results in a slower decrease in the oxygen abundance.
We have assumed ℓ =   ≈ 10 8 cm in all our calculations.While this is likely a reasonable assumption for overturning convection, in the thermohaline regime the mixing length is not well-constrained.We have shown analytically in section 2.4 that, for the limit Γ << 1 we can reproduce the thermohaline mixing parameter  =   ℓ derived in previous studies (e.g.Kippenhahn et al. 1980).Montgomery & Dunlap (2024) used the thermohaline prescription of Brown et al. (2013) to obtain where  is the kinematic viscosity and  is a constant.They also used a different scaling  ∝  −0.62 0 which they found gave a better fit to the numerical data.
These different prescriptions for  are compared in Figure 6 using the solutions found for Γ, as well as the relevant gradients for the input parameters given by our WD model of 0.9  ⊙ .To study the dependency on ℓ we considered two choices of mixing length: ℓ =   , as in the overturning regime, and ℓ = 100 cm as suggested by Montgomery & Dunlap (2024).Interestingly, we find that the values of  agree within an order of magnitude with the value directly computed from the solution of Γ as  =   ℓ = 2 0   Γ.This may be related to the fact that, in the thermohaline regime, the value of the product   ℓ ∝ Pe and therefore the mixing coefficient is set directly by the composition flux as shown in F23.Although different prescriptions have different predictions for   and ℓ individually, their product is determined once the composition flux is specified.
As a consequence of the almost constant thermohaline mixing coefficient, the convective velocity is affected by the mixing length assumption as roughly   ∝ ℓ −1 , which is why F23 found such a strong enhancement of   by rotation, which reduces the effective mixing length.For ℓ = 100 cm, the convective velocity in the thermohaline regime is increased by up to six orders of magnitude compared to what is shown in Figure 4. Despite this large increase,   remains ≲ 0.5 cm s −1 for most of the evolution, since we previously found that   has typical values ≲ 10 −6 cm s −1 ( ≲ 10 −1 after a few Myr) when assuming ℓ =   .

Summary and Conclusions
We have developed a general mixing length theory of convection (MLT) that self-consistently includes the effects of thermal diffusion and composition gradients.Our formulation smoothly transitions between the regimes of fast adiabatic convection at large Peclet number and slow thermohaline convection at low Peclet number.It also allows for both thermally driven and compositionally driven convection, including correctly accounting for the direction of heat transport for compositionally driven convection in a thermally stable background.Once implemented into 1D evolution codes, this will enable investigations of the compositionally driven mixing that occurs in white dwarfs during crystallization.
Our results improve on existing implementations of MLT in stellar evolution codes.Rather than including the effects of composition by adding composition gradients into the existing MLT equations (with different codes taking different approaches, see discussion in Section 2.2), we instead use the heat and composition fluxes as input quantities to derive the resulting convective velocity and thermal and compositional gradients.This leads to a 5th order polynomial (Equation ( 22)) that can be solved for the convective efficiency or convective velocity, and from there for the thermal and composition gradients (Equations ( 11), (12), and ( 21)).This expression replaces the usual third-order polynomial used in stellar evolution models (Equation (10)).A similar set of equations was found by Umezu & Nakakita (1988); Nakakita & Umezu (1994); Umezu (2008Umezu ( , 2009) ) in the context of core convection in massive stars.
In stellar evolution codes, our MLT description could be implemented using the heat and composition fluxes obtained from the conservation laws for energy and matter.For example, in the case of white dwarf crystallization, the energy flux depends on the internal energy of the mixture, the latent heat, and the gravitational energy released.The composition flux depends on the correct modeling of the chemical separation during the growth of the crystallized core.Once the composition flux is determined, our MLT prescription can be used to find the corresponding composition gradient for the next time step, similar to the way in which the heat flux is used to determine the temperature gradient.
In this study, we did not attempt to implement our theory in an existing or new code.However, as a first step in applying this theory, we use the output of fiducial WD cooling models from MESA to investigate the regimes of convection that would appear in crystallizing WDs.In doing so, we calculated the heat and composition fluxes above the crystallization front at every stage of the cooling and used them to compute the properties of convection during the star's evolution.We found that, during most of the cooling evolution, convection is in the thermohaline regime, characterized by a small efficiency and convective velocities   ∼ 10 −9 -10 −5 cm s −1 .However, at the onset of the crystallization, all the models are found in the efficient regime for a short time (≲ 10 Myr), with much larger values of   ∼ 10-100 cm s −1 .This happens because the growth rate of the solid core is large at early times, and its small size also leads to large values of the composition flux (exceeding the dimensionless composition flux  = 1, Equation ( 21), which leads to fast convection).
These results have interesting implications for white dwarf magnetic fields.The efficiency of compositionally driven convection and the corresponding convective velocities play an important role in the crystallization-driven dynamo model proposed to explain highly magnetized WDs (Isern et al. 2017;Ginzburg et al. 2022).For a saturated dynamo with  2 /4 ∼  2  (e.g.Isern et al. 2017), the magnetic field spans the range  ∼ 10 −6 -10 5  1/2 6 G for the range of convective velocities that we find above,   ∼ 10 −9 -10 2 cm s −1 (assuming a typical density of 10 6 g cm −3 ).It is important to note that our estimates of the convective velocity do not take into account rotation or magnetic forces.As shown in numerical simulations and scaling laws from mixing length theory, those effects can significantly change the magnitude of the convective velocity (e.g., Mochkovitch 1983;Ginzburg et al. 2022;Fuentes et al. 2023).However, it is intriguing that the convective velocities we find in the very early phases of crystallization appear to be large enough to generate magnetic fields comparable to those observed.We explore the possibility of a short-lived intense dynamo following the onset of crystallization and the expected magnetic field strengths taking into account rotation in a companion paper (Fuentes et al. 2024) Our calculations can be improved in several respects.The white dwarf cooling models that we use include chemical separation following Bauer (2023), in which it is assumed that the fluid layer ahead of the crystallization front is fully mixed, and the gravitational energy release from chemical separation is deposited by hand in the model.Implementing the MLT derived here in a time-dependent model would allow to follow the heat and compositional transport in detail, and solve the relevant gradients to construct more consistent tempera-ture and composition profiles during the cooling.The onset of crystallization and the short-lived phase of fast convection needs to be followed in more detail, and its dependence on white dwarf mass and the extent of the convection zone during this time determined (which depends on the C/O profile and sets the delay before the field emerges at the white dwarf surface; Blatman & Ginzburg Also, as mentioned above, rotation and magnetic forces should be included when deriving the convective velocity.F23 estimated that rotation will increase the convective velocity significantly in the thermohaline regime (although the convection is still inefficient, with low Peclet number) (see also Mochkovitch 1983;Ginzburg et al. 2022).
More sophisticated transport models beyond MLT are needed, in particular in the thermohaline regime.For example, Brown et al. (2013) derived expressions for the fluxes in thermohaline convection based on a series of numerical simulations covering a range of values of the governing parameters.A recent paper by Montgomery & Dunlap (2024) used these results (with some modifications) to study mixing in white dwarfs following crystallization.In agreement with F23, they concluded that the convective velocities are much lower than estimated by Isern et al. (1997) and Ginzburg et al. (2022).However, they also advocated for a much smaller mixing length (ℓ ∼ 100 cm) than we have assumed here (ℓ ≈   ∼ 10 8 cm).To investigate this point, we compare the diffusivities from Brown et al. (2013) with the MLT results for both of these choices of ℓ in Section 4. We find that the mixing coefficients are within an order of magnitude regardless of the chosen prescription or mixing length.The reason for this is that, in the thermohaline regime, the Peclet number ∝   ℓ is set by the composition flux ( ≪ 1) (as shown by F23), therefore is the effective mixing coefficient which also scales as   ℓ.However, given this result, the convective velocity increases for a smaller ℓ, scaling as ℓ −1 .This is consistent with the enhancement of   found by F23 in the rotating case as the effective mixing length is reduced.
Further multidimensional numerical simulations are also needed to study the transport properties of thermohaline convection under conditions of white dwarf interiors.For example, recent simulations have shown that background magnetic fields can significantly enhance the rate of chemical mixing by fingering convection (e.g., Harrington & Garaud 2019;Fraser et al. 2024).In this context, it would be also interesting to know how fingers in the thermohaline regime react back on a large intensity field created at the onset of crystallization.
We thank Adrian Fraser for useful conversations and pointing toward important references regarding thermohaline convection. A. -10 3 -10 0 -10 3 0 10 3 10 0 10 3 -10 2 -10 0 -10 2 -10 4 B. The Other Solution for  < 0 and ∇ rad /∇ ad > 0 As discussed in Section 3, the quintic Equation ( 22) has a second solution in the region of parameter space  < 0 and ∇ rad /∇ ad > 0. For completeness, we show the second solution in Figure 8, which should be compared with Figure 1.Nakakita & Umezu (1994) showed that the second solution is unstable for some choices of mixing length, but were not able to completely rule out that it could be relevant for overturning convection in massive stars (see also Umezu 2009).Here, looking broadly across the parameter space, we see that Γ changes discontinuously across  = 0 in the overturning convection regime.Since we do not expect the convective velocity to change drastically on going from a small inward composition flux to outward composition flux, we reject this solution as unphysical.In addition, this solution has the property that Γ decreases with increasing composition or heat flux, which is opposite to the expected behavior.

Figure 1 .
Figure 1.The convective efficiency Γ (Equation (3)) as a function of the total heat flux as measured by ∇ rad /∇ ad (Equation.(8)) and the composition flux as measured by  (Equation (21)).The solutions were obtained by solving Equation (22) numerically.We show the different regimes of convection, as well as the stable region where convection does not occur.The red line shows the transition to the thermohaline convection regime and the magenta line delimits the regime of semiconvection.To solve Equation (22), we used  0 = 9/4,  = 2.5 × 10 4 , and ∇ ad = 1/3.

Figure 2 .
Figure 2. Normalized temperature and composition gradients ∇/∇ ad (left panel) and ∇ com /∇ ad (right panel) as a function of ∇ rad /∇ ad and .The solutions were obtained using equations (11) and (21), and the Γ values from Figure1.In regions that are stable to convection, we set Γ = 0 and ∇ = ∇ rad .

Figure 3 .
Figure 3. Convective heat flux   (normalized by  ad =      ∇ ad /  ) as a function of ∇ rad /∇ ad and .These values were obtained using Equations (18) and (19), and the Γ values from Figure 1.

Figure 4 .
Figure 4. Evolution tracks in the -∇ rad /∇ ad parameter space for three white dwarf cooling models.The solid lines indicate the models with initial abundances from stellar evolution, while the dashed line denotes the model with a uniform 50/50 C/O core abundance.In this diagram, the star moves from right to left as it evolves.The background color map indicates the magnitude of the convective velocity   computed from Γ.We use  0 = 9/4,  = 10 11 , ∇ ad = 1/3,   = 50 cm 2 s −1 , and ℓ =   = 10 8 cm.

Figure 5 .
Figure 5.Time evolution of the convective parameters above the solid-liquid interface following the onset of crystallization for the models shown in Figure 4. To follow the crystallization front and the flux of light elements out of the solid core, we used Equation (30), the output of the Skye EOS of Jermyn et al. (2021), the phase diagram of Blouin et al. (2020), and the phase separation routine in Bauer (2023) (see text for details).

Figure 6 .
Figure6.A comparison of different prescriptions for the thermohaline mixing coefficient.We used the stellar parameters and values of Γ, ∇, and ∇ com obtained with our MLT for the WD model of 0.9  ⊙ to compute the different scaling relations with  0 .We assumed  0 = 9/4 for the relation shown in Equation (26) (blue), and  = 2  and  = 0.4 cm 2 s −1 for the relations found byMontgomery & Dunlap (2024) (orange and green).Additionally, we compare the results with the expected diffusion parameter from  =   ℓ = 2 0   Γ for two different values of  0 (red and magenta).
C. acknowledges support by NSERC Discovery Grant RGPIN-2023-03620.J.R.F. is supported by NASA through grants 80NSSC19K0267 and 80NSSC20K0193.A.C. M.C.-T.are members of the Centre de Recherche en Astrophysique du Québec (CRAQ) and the Institut Trottier de recherche sur les exoplanètes (iREx).
, .The subscript  is used to represent the partial derivatives taken at constant composition.The temperature gradient with pressure in the star is ∇ =  ln / ln | ★ and the temperature gradient with pressure experienced by a fluid element as it moves is ∇  .The buoyancy force associated with the density excess   gives the convective velocity as1  2  ≈ −( /) (ℓ/4) or :   and  are the density and temperature excess carried by a fluid element, respectively, ℓ is the mixing length,   is the pressure scale height,   =  ln / ln  | , , and   =  ln / ln |