Fluid Mixing during Phase Separation in Crystallizing White Dwarfs

Accurate models of cooling white dwarfs must treat the energy released as their cores crystallize. This phase transition slows the cooling by releasing latent heat and also gravitational energy, which results from phase separation: liquid C is released from the solid C/O core, driving an outward carbon flux. The Gaia color–magnitude diagram provides striking confirmation of this theory by revealing a mass-dependent overdensity of white dwarfs, indicating slowed cooling at the expected location. However, the observed overdensity is enhanced relative to the models. Additionally, it is associated with increased magnetism, suggesting a link between crystallization and magnetic field generation. Recent works aimed at explaining an enhanced cooling delay and magnetic field generation employ a uniform mixing prescription that assumes large-scale turbulent motions; we show here that these calculations are not self-consistent. We also show that thermohaline mixing is most likely efficient enough to provide the required chemical redistribution during C/O phase separation, and that the resulting velocities and mixing lengths are much smaller than previous estimates. These reduced fluid motions cannot generate measurable magnetic fields, suggesting any link with crystallization needs to invoke a separate mechanism. Finally, this mixing alters the chemical profiles, which in turn affects the frequencies of the pulsation modes.


ASTROPHYSICAL CONTEXT
The Gaia mission has uncovered a clear sequence in the white dwarf color-magnitude diagram that confirms theoretical expectations (van Horn 1968) that the steady cooling of white dwarfs (WDs) slows as they crystallize (Tremblay et al. 2019).Termed the "Q-branch," this pile-up of high mass (0.9 -1.1 M ⊙ ) WDs coincides with carbon/oxygen models that have begun crystallizing.However, the detailed shape of the observed bump in the WD luminosity function is both narrower and higher than predicted by the models, indicating an extra cooling delay, possibly up to 8 Gyr for some fraction of the population (Cheng et al. 2019).This has led to renewed interest in the crystallization of mixtures and in 22 Ne sedimentation and "distillation" (Blouin et al. 2020(Blouin et al. , 2021;;Bauer et al. 2020).
Numerous calculations, using molecular dynamics (Horowitz et al. 2010;Schneider et al. 2012), density functional theory (Segretain & Chabrier 1993), and free energy approaches (Blouin et al. 2020;Medin & Cumming 2010), support the claim that the different ionic species in the interior of a white dwarf should partially separate upon crystallization.For a binary elemental mixture, the high-Z species should be enhanced in the solid relative to the low-Z species, with the surrounding fluid being depleted in the high-Z element.Subsequently, through a Rayleigh-Taylor (RT) instability, the lower-Z (and therefore lighter) liquid layer is assumed to mix with the denser layers above it (e.g., Isern et al. 1997;Salaris et al. 1997).This chemical redistribution releases gravitational energy, providing an additional energy source which slows the cooling of the WD (e.g., Stevenson 1980;Mochkovitch 1983;Isern et al. 1997;Salaris et al. 1997;Montgomery et al. 1999;Althaus et al. 2003;Bauer 2023).This mixing also systematically alters the g-mode oscillation frequencies of these stars (e.g., Montgomery & Winget 1999;Córsico et al. 2005;De Gerónimo et al. 2019).

MIXING RE-EXAMINED
The suggestion has recently been made that the fluid motions resulting from phase separation could be large enough to generate magnetic fields of order ∼ 0.1 MG (Isern et al. 2017), and even larger fields are estimated in the recent calculations of Ginzburg et al. (2022).We examine the viability of this hypothesis through estimates of the mass flux and fluid velocity required in the standard treatment of phase separation.We examine different mixing mechanisms and discuss the consequences our results have for magnetic field generation and pulsations in these stars.
From the perspective of a 1-D stellar evolution code, at the end of a given time step the model will have increased its crystallized mass fraction, leaving a thin carbon-enhanced layer on top of the crystallized core.Since carbon is lighter than oxygen at these conditions ((δρ/ρ) 0 ∼ 3 × 10 −3 ), this fluid layer will be RT unstable to mixing.Given the large acceleration due to buoyancy felt by a pure carbon fluid el- The evolution of the oxygen profile due to crystallization and phase separation for a 1 M⊙ C/O WD.The different curves are labeled with the T eff , percent crystallization by mass, and age of the model, respectively.Initially, the C/O core has a uniform composition out to Mr ∼ 0.9 M⋆, with an oxygen mass fraction of 0.5.The dotted lines indicate the boundary between the solid and liquid regions in each model.The carbon profile for Mr/M⋆ < 0.9 is essentially equal to 1 − XO. (b): The carbon flux (blue solid curve, left axis) and velocity times carbon enhancement (red dashed curve, right axis) of carbon-enriched bubbles derived from two consecutive evolutionary models at ∼ 16 % crystallization by mass.ement (g δρ/ρ ∼ 10 5 cm/s 2 ), this mixing is assumed to be rapid and vigorous, producing a layer of uniform composition out to a point in the model where the carbon mass fraction matches that of the layer above it.This ensures that the resulting configuration is stable to further RT mixing.In section 5 we explore the case in which the carbon enhancement of a fluid bubble is much smaller, resulting in correspondingly smaller accelerations.
The detailed mechanics of this mixing are not usually considered, e.g., the size or velocity of the convective bubbles.Recently, Isern et al. (2017) have examined this mixing in more detail.They assumed bubbles of radius ∼ 0.1 R core , where R core is the radius of the crystallized core.Balancing the buoyancy force with that due to turbulent drag, they computed velocities of ∼ 35 km s −1 .Using relations in Christensen (2010), they suggested that this could lead to a dynamo process that could generate magnetic fields of order 0.1 MG.More recently, Ginzburg et al. (2022) have estimated a much smaller carbon enhancement for the fluid bubbles, resulting in reduced mixing velocities of order 10 2 cm s −1 .
From evolutionary calculations we know that the rate of crystallization should be small, with WDs taking ∼ 10 9 years to crystallize most of the mass in their cores; at any given point in the crystallization process, the radius of the crystallized core is increasing at a rate of approximately 10 −8 cm s −1 .Therefore, the thickness of the carbonenhanced layer above the crystallized core should be increasing at a similar rate.Thus, we expect that the net outward flux of carbon should be very small.We estimate these numbers in the following section.

CHEMICAL TRANSPORT IN THE STANDARD PICTURE
For the calculations described in this and subsequent sections we have used version 15140 of the stellar evolution code MESA.In particular, the numerical details of crystallization and fluid mixing, both neutrally buoyant and the standard prescription, are described in Appendix A. We note that more recent versions of MESA have the ability to do the standard uniform mixing during crystallization as a built-in option (Bauer 2023), but we have not used this functionality.
We calculate results of the standard prescription for crystallization and phase separation in an evolving WD model, assuming that the fluid mixing leads to a uniformly mixed overlying fluid layer.This model has a carbon/oxygen core and a total mass of 1 M ⊙ .We have artificially chosen an oxygen profile that extends past 0.9 M ⋆ in an effort to maximize the fluid mixing during subsequent evolutionary stages and a high mass so that the model will be partially crystallized in the DAV instability strip.For the phase diagram of carbon/oxygen crystallization, we use the recent results of Blouin et al. (2020).In Figure 1a we show the evolution of the chemical profile of oxygen as the model cools.The apparent discontinuity in the profiles (dashed portion of curves) indicates the edge of the crystallized core and the fluid layers above it.
Using consecutive models with differing crystallized mass fractions, we can compute the mass flux of carbon as a function of radius that is necessary to transform the carbon profile at one time step into the carbon profile at the next time step; we plot this as the solid blue curve (left axis) in Figure 1b.We note that when the model is about 14% crystallized by mass, the uniformly mixed fluid region extends from ∼ 0.1 to ∼ 0.9 in M r /M ⋆ .This stage of crystallization was chosen in an effort to maximize the extent of the mixed region and therefore the fluid velocities.
Next, we compute the magnitude of fluid velocities that are required to produce this mixing.Assuming that the carbon enhancement of a fluid bubble relative to its environment is given by ∆X C , and that the bubbles have a filling fraction of Figure 2. The same as Figure 1, except for the case of neutrally buoyant mixing.(a): Instead of flat profiles in the mixed regions, these profiles have a non-zero slope that results in N 2 ∼ 0. (b): Using equation 9 we are able to solve for the velocity times the mixing length, which is shown on the right-hand axis, again assuming f = 0.5.We note that vC × L ∼ 1 cm 2 / s.
f , then the flux of carbon, F C , is related to the velocity, v C , by where ρ is the mass density.Taking f = 0.5, we find that the product v C × ∆X C is given by the red dashed curve (right axis) in Figure 1b.If ∆X C ≈ 10 −2 , then the velocities are quite low, of order ∼ 10 −8 cm s −1 , which is ∼ 10 14 times smaller than the estimate of Isern et al. (2017).On the other hand, if we assume values of ∆X C of 10 −16 and 10 −12 , then we obtain velocities similar to those of Isern et al. (2017) and Ginzburg et al. (2022), respectively.Thus, the velocities needed to generate appreciable magnetic fields or achieve turbulent mixing require very small values for ∆X C .However, in the following section, we show that such a small carbon enhancement results in bubbles that are limited to small vertical displacements.
4. PHYSICAL INCONSISTENCY The mixing envisioned in Ginzburg et al. (2022) and Isern et al. (2017) is vigorous enough to produce a uniform chemical profile in the mixed region and is relatively rapid in the sense that the convective turnover time, τ c , is short compared to a bubble's thermal diffusion time scale, τ * KH (see Kippenhahn et al. 1980).To see this, we use where L is the radius of a spherical bubble and κ T is the thermal diffusivity as defined in Appendix C. Through most of its evolution, the assumed mixing length in the model of Ginzburg et al. is given by a pressure scale height.Taking L as half of this mixing length, we find a thermal diffusion time scale for the bubble of τ * KH ∼ 2 × 10 6 yr.The ratio of τ c (∼ 0.1 yr) to τ * KH is, then, ∼ 10 −7 , meaning that there will be almost no heat exchange between the bubble and its environment.Thus, the thermodynamics of the bubble will essentially be adiabatic; yet, as we will show, this contradicts the assumptions made in deriving τ c .
In uniform composition convection zones in stellar atmospheres, the temperature decreases outwards so that when a fluid parcel rises and expands in response to the lower pressure, the possibility exists for it to be less dense than its surroundings and continue to rise.This will happen if, when it cools adiabatically, it is still hotter than its surroundings (this is just the Schwarzschild criterion).In contrast, in a nearly isothermal environment (still with uniform composition) such as the core of a WD, a rising element that expands adiabatically will be cooler than its surroundings, making it denser and hence convectively stable.On the other hand, if the composition of the rising parcel is such that it is less dense than the surrounding fluid, then it may continue to rise, but the cooling effect on its density must still be taken into account.The derivation of convective velocities in Ginzburg et al. (2022) only considers the effect of composition on density and not the effect of adiabatic cooling on density, which can very quickly dominate in a nearly isothermal environment.We quantify this argument below.
The return force on an adiabatically displaced fluid element in a uniform composition medium is proportional to N 2 δr, where δr is the vertical displacement of the fluid element from its equilibrium position and N 2 is the square of the Brunt-Väisälä frequency.On the other hand, if the fluid element is enhanced in carbon relative to its environment, it will experience an additional upward buoyancy force proportional to g(δρ/ρ) 0 ∆X C , where (δρ/ρ) 0 is the fractional density contrast between a pure oxygen and a pure carbon plasma (∼ 3 × 10 −3 in our models), and ∆X C is the enhancement of the mass fraction of carbon in the fluid element compared to its surroundings.Including both effects yields the following differential equation for the acceleration of the fluid element: Since the background chemical profile is assumed to be uniform, ∆X C is independent of the height δr, so the second term on the RHS of equation 3 is approximately constant.On the other hand, the first term is proportional to δr, so it becomes more negative the farther up the fluid element travels.Ginzburg et al. (2022) omit the first term on the RHS of equation 3 in their analysis, so their fluid elements can rise without bound.This naturally leads them to assume a mixing length that is of order a pressure scale height.However, when we include both terms in equation 3, we find that a fluid element that starts out at δr = 0 with a positive carbon enhancement will eventually rise to a new equilibrium point where these two terms balance each other.This equilibrium distance is given by In other words, the value of δr eq in equation 4 is the actual mixing length L that should be used based on these assumptions.
Using the values from Ginzburg et al. ( 2022) we estimate their model predicts (δρ/ρ) 0 ∆X C ∼ 10 −14 , and for values of g and N 2 , we use a fiducial WD model.We find that L ∼ 0.03 cm.This value is approximately 10 orders of magnitude smaller than their assumed mixing length of ∼ 10 8 cm, highlighting a severe inconsistency.
To summarize this chain of reasoning, the assumptions of Ginzburg et al. (2022) of (1) a large mixing length (of order a pressure scale height), (2) a bubble size that is of order this mixing length, and (3) instantaneous thermal equilibration of the rising fluid element with its surroundings (i.e., neglect of the adiabatic term in equation 3), lead to fluid velocities that imply a convective turnover time that is much shorter than the thermal diffusion time scale.This requires the resulting fluid motions to be adiabatic, which is contrary to (3).Thus, the model of Ginzburg et al. (2022) is invalid in this context.

NEUTRALLY BUOYANT MIXING
In this section, we derive conditions whereby weak RT mixing could produce a profile that is very close to neutral stability.By assuming that the fluid motions are adiabatic, we are implicitly ignoring the physics related to thermohaline mixing, which we examine in detail in sections 7 and 8.The main motivation for this section is to provide a calculation of an upper limit for the steepness of the chemical profile.
The Brunt-Väisälä frequency can be written in a form that explicitly displays its dependence on the composition profile: where we have written N T and N µ for the thermal and compositional parts of the Brunt-Väisälä frequency, respectively.In the above formulae, P is the pressure, where ⃗ X k is the vector of mass fractions of the different chemical species defined at the k th mesh point.A neutrally buoyant profile (N 2 ≃ 0) will be one for which B = ∇ − ∇ ad ; this specifies the required gradient in the carbon (and oxygen) profile.
We show the result of neutrally buoyant mixing in Figure 2a.For the sequence of five models shown, the mixed region now has a positive oxygen gradient such that N 2 = 0. Also of note, unlike the uniform mixing case in which the mixing region extends out to M r /M ⋆ ≈ 0.9 even for small amounts of crystallization, here the mixed region is smaller, only reaching the edge of the core as the crystallized mass fraction itself moves outward.
In Figure 2b, we show the carbon flux required to produce this flux.Independent of any particular mixing length theory, the carbon enhancement of a fluid bubble, δX C , will be the difference in composition between the region in which it forms and the region in which it eventually mixes.Thus, if the fluid element travels a distance L, this enhancement is related to the C profile by yielding, with equation 1, the following equation for the flux: (for an essentially identical result, see equation 19 of Mochkovitch 1983).Within a factor of order unity, the Reynolds number, Re, of this flow is approximately given by Re = v C L/ν kin , where the kinematic viscosity, ν kin , is approximately 0.4 cm 2 s −1 for these conditions (Nandkumar & Pethick 1984).In Figure 3, we show the value of Re for this model.We see that Re is of order unity, which in general is much too small for the flow to be turbulent; typically, Re > 10 3 is required for turbulent flow. 1 While the idea of "laminar mixing" may sound like a contradiction of terms, we note that laminar motions have been considered as part of the process for mixing heterogeneities in the Earth's mantle (e.g., Olson et al. 1984).
6. PULSATIONS AND NEUTRAL PROFILES Montgomery & Winget (1999) were the first to treat in detail the effect that a crystallized core would have on the g-mode pulsations of WD stars.They found that g-modes are excluded from the crystallized solid region of the models.Due to this shrinking of the g-mode cavity, the average period spacing of consecutive radial overtones of a given ℓ, ⟨∆P ⟩ ℓ , was found to be a monotonically increasing function of the crystallized mass fraction.For instance, in going from uncrystallized to 90% crystallized by mass, while otherwise holding the thermal and mechanical structure constant, ⟨∆P ⟩ ℓ=1 increases from 31 s to 43 s.
Since g modes locally propagate in regions in which their squared frequency is less than N 2 , they will be excluded from neutrally mixed regions with N 2 = 0 as well as the crystallized core.In this case, we expect the effect of the mixed plus crystallized region to mimic that of a star with a larger crystallized core.As we see in Figure 4, when our fiducial model is 20% crystallized by mass it will have a crystallized plus mixed region that is almost 50% of the total mass.Thus, if we fit the pulsation frequencies of this model with a model that does not take into account neutrally buoyant mixing, the fit will erroneously favor models that are 50% crystallized by mass.In fact, one prediction of this model is that stars analyzed with traditional seismic fits will be found to be either 0% crystallized or ≳ 40% crystallized by mass because the rapid growth of the neutrally buoyant region will make it unlikely to find them inside of these limits.These general considerations should apply to all models which are massive enough to be partially crystallized while in the instability strip (e.g., Montgomery & Winget 1999;Córsico et al. 2005;De Gerónimo et al. 2019).
In the following section, we examine the potential effects of thermohaline mixing on the chemical profiles.If the thermohaline diffusion coefficient, D thrm , is sufficiently large, the resulting profiles will not be as steep as in the neutrally buoyant case.The associated pulsational models will, therefore, also be different from those considered here.However, the simulations used to derive D thrm (e.g., Brown et al. 2013) must be extrapolated by ∼ 5 orders of magnitude in τ ≡ κ µ /κ T (the ratio of molecular to thermal diffusivity) and 3 to 4 orders of magnitude in ν kin /κ µ (i.e., the Schmidt number) to match the WD interior conditions considered here (see Table 1 for definitions and values of these variables).Given this mismatch, neutral profiles cannot be definitively ruled out.In addition, asteroseismology could, in principle, provide evidence for the efficiency of the thermohaline mixing discussed below.

THERMOHALINE MIXING
The treatment of mixing in section 5 assumes that the fluid elements do not exchange heat with their environment, i.e., the motion is adiabatic.In this case, fluid motions occur only when the chemical gradients are large enough that N 2 < 0. However, when heat exchange with the environment is considered, smaller chemical gradients can produce mixing in a process termed thermohaline convection or thermohaline mixing.If this is sufficiently efficient, it could provide the required mixing during phase separation and prevent the formation of steep chemical gradients that would produce an RT instability.
In this section, we use MESA's thermohaline capabilities to see if the predicted fluxes of C and O due to thermohaline mixing are large enough to fulfill the mixing requirements of our models.In order to perform this calculation, we have modified MESA so that its thermohaline prescription uses values for the kinematic viscosity ν kin (Nandkumar & Pethick 1984) and the molecular diffusivity of the ions κ µ (Caplan & Freeman 2021;Caplan et al. 2022) that are appropriate for WD stars.For these models we find κ µ ≈ 10 −5 cm 2 s −1 and ν kin ≈ 0.4 cm 2 s −1 (for these and other values, see Table 1 in Appendix C).
To isolate this mixing from convective mixing, which would arise if at any point during convergence the profiles were to become too steep, we construct a chemical profile that is only 90% as steep as a neutrally buoyant profile, i.e., a profile with B = 0.9 (∇ − ∇ ad ), since we can study this state in the absence of convection.Then, in the following time step, we turn on thermohaline mixing and observe the flux of carbon that this mixing produces.At this stage we have to be careful to reduce the time step so that the profile is not appreciably altered during the time step.Since the flux scales as the chemical gradient to the −1.62 power (R −1.62 0 , see equation B14 in section B.1), we divide the flux by 0.9 1.62 ≈ 0.84 to estimate the flux that a neutrally buoyant profile would produce.The results of this procedure are plotted in Figure 5.The blue curve shows the flux of C that is required to maintain this profile as computed from two consecutive time steps.The orange curve shows the net flux of C due to thermohaline mixing that this profile would produce when using the formalism of Kippenhahn et al. (1980), while the green and red curves show the predicted fluxes when using the Brown et al. (2013) and Traxler et al. (2011) prescriptions, respectively.Since the fluxes due to the Kippenhahn et al. (1980) and Brown et al. (2013) prescriptions exceed that needed to produce and maintain a neutrally buoyant profile, for these theories thermohaline mixing would be efficient enough to produce all the required mixing and the RT instability would not come in to play.On the other hand, the Traxler et al. (2011) prescription produces mixing which is too weak to generate the required fluxes; in this case, we suggest that a neutrally buoyant profile would form and the RT instability described in the section 5 would operate.

A NEW MIXING LENGTH THEORY FOR THERMOHALINE MIXING
Given that the Brown et al. (2013) prescription is the most recent, is derived from an underlying physical model, and has been benchmarked against numerical simulations, we treat it as the currently favored theory for thermohaline processes.From Figure 5, we see that it can transport a factor of 40 or more carbon than the models require, making it the dominant process for chemical transport during C/O phase separation.While we have not self-consistently computed the resulting chemical gradients, we can estimate that they will be reduced by a factor of 40 1/1.62 ≈ 10 from the neutrally buoyant case.Thus, the profiles will be closer to the flat, uniformly mixed profiles of the standard case than the neutrally buoyant case.
On the other hand, Brown et al. (2013) find somewhat less good agreement between their theory and the simulations for the astrophysically interesting low-Pr, low-τ regime.This led us to re-examine their model, which we present in Appendix B.1.Then, in Appendix B.2, we show that its general form can be derived as a mixing length theory and we suggest a small addition that accounts for the loss of compositional flux due to turbulent diffusivity.This decrease in efficiency of the chemical transport allows us to provide a closer match to their low-diffusivity, low-viscosity numerical simulations (see Figures 6a,b).
Our new predicted fluxes are about a factor of two smaller than those of Brown et al. (2013), and this is not large enough to qualitatively change our conclusions in section 7. We estimate that the resulting chemical profiles should produce a 16% reduction in N 2 compared to the uniform mixing case, a signature with potential asteroseismological consequences.

MAGNETIC FIELD GENERATION
The high-mass overdensity defining the Q-branch contains an enhancement in the incidence of magnetic white dwarfs (McCleery et al. 2020).Among lower-mass WDs (M ⋆ ≲ 0.75 M ⊙ ), Bagnulo & Landstreet (2021, 2022) find an increasing incidence of magnetic fields as they cool, particularly after the point when models of these stars predict the onset of crystallization.In addition, Parsons et al. (2021) note a lack of identifiable progenitor systems for magnetic cataclysmic variables.This leads Schreiber et al. (2021Schreiber et al. ( , 2022) ) to speculate that the phenomenon of crystallization, which occurs at late times, may be crucial for the generation of these fields.These groups cite Isern et al. (2017) and/or Ginzburg et al. (2022) as a possible explanation for this.In a novel use of Isern et al. (2017), Camisassa et al. (2022) infer the existence of O/Ne cores based on the existence of magnetic fields in white dwarfs too hot to be crystallizing if they have C/O cores, but which are expected to be crystallizing if they are O/Ne.
Based on work by Augustson et al. (2016Augustson et al. ( , 2019)), Ginzburg et al. (2022) use the following equation to estimate the strength of the dynamo-generated magnetic fields, where P is the rotation period of the model.In their models they find v c L ∼ 10 10 cm 2 s −1 .However, we show in section 4 that their approach is not physically self-consistent, i.e., their assumption of rapid, adiabatic mixing is not consistent with an assumption of large mixing lengths.In our neutrally buoyant mixing prescription, the carbon enhancement of the bubbles is tied to the carbon gradient in the neutrally buoyant regions, so it is not possible in our model to have both large mixing lengths and large velocities.Instead, we find v c L ∼ 1 cm 2 s −1 .For the thermohaline case, we find values that are slightly larger, with v c L ∼ 200 cm 2 s −1 .Scaling from the results of Ginzburg et al. (2022) to our thermohaline results, we predict magnetic fields ∼ 10 4 times weaker, i.e., B ≲ 10 kG.We conclude that the fluid mixing associated with phase separation is not a viable mechanism for generating large magnetic fields in WDs.However, this does not mean that crystallization itself cannot still play a role in magnetic field generation.If a WD is differentially rotating, either due to prior evolution or active accretion of angular momentum, then its rotational profile will be altered by crystallization.The crystallized core will rotate as a solid body while the overlying fluid layers will not.Thus, a strong shear layer could be set up at the solid/fluid boundary, and the fluid motions associated with this region plus the rotation of the WD could lead to magnetic field generation (e.g, Spruit 2002).
10. DISCUSSION Given the long time scale for crystallization, we have reconsidered the assumption of other studies that the mixing due to phase separation is vigorous and turbulent.If thermohaline mixing is too inefficient, we assume that slow RT mixing will produce neutrally buoyant profiles having N 2 = 0, a result found in recent studies (Brooks et al. 2017;Schwab & Garaud 2019), rather than the uniformly mixed profiles that are normally assumed.Simple estimates show that v c L ∼ 1 cm 2 s −1 , with a Reynolds number of O(1).
On the other hand, current theories of thermohaline convection predict that regions with shallower profiles than neutrally buoyant are capable of sustaining thermohaline fluxes large enough to transport the carbon released during phase separation.Assuming the Brown et al. (2013) formalism extrapolates to the physical regime considered here, we find that v c L ∼ 200 cm 2 s −1 and a Reynolds number of ∼ 500.The theory also allows us to independently estimate the mixing length as the approximate length of the thermohaline "fingers" (∼ 100 cm) and the velocities as ∼ 2 cm s −1 .
In either the neutrally buoyant or thermohaline mixing case, we find that these mixing lengths and velocities are much smaller than the estimates in Ginzburg et al. (2022) and Isern et al. (2017).In fact, the velocities and vigorous turbulent convection of those models result in physical inconsistencies; we thus conclude that such mixing is not a viable mechanism for generating large magnetic fields in WDs.
While we have confined our discussion to WDs with carbon/oxygen cores, everything we have presented should apply equally well to WDs with crystallizing O/Ne/Mg cores.First, for the reasons cited in the above paragraphs, the fluid motions will not be vigorous enough to produce measurable magnetic fields in these WDs.Second, for the neutrally buoyant mixing case, the mapping of M cryst in the standard treatment to M cryst + M mix (as shown in Figure 4) should still hold.Although the cores of these models may be mostly crystallized by the time the stars reach the DAV instability strip, for high-mass DBs our results could be important for interpreting their pulsations as they begin to crystallize in the DBV instability strip.
Finally, we anticipate that implementing a thermohaline mixing prescription may be relevant for other cases involving phase separation and mixing, such as the process of 22 Ne distillation in ultramassive WDs as proposed by Blouin et al. (2021).

APPENDIX A. NUMERICAL IMPLEMENTATION OF NEUTRALLY-BUOYANT MIXING
We implement this new mixing prescription in the stellar evolution code MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019)).MESA provides customization for different physical processes through the run_star_extras.ffile, which allows for a large range of user-specified modifications without the need to modify any central routines.The mixing is applied at the end of each time step in the function extras_finish_step, and heating due to the release of gravitational energy is then applied during the next time step through the other_energy hook in extras_controls.These calculations were performed using versions 12778 and 15140 of MESA.
We calculate N 2 and B as given in equations ( 5) and ( 7) by calling do_brunt_N2 and do_brunt_B in MESA's brunt module.We adjust the size of the mixed region based on the amount of carbon that is liberated from the newly crystallized shells.While our procedure does not result in an N 2 which is identically zero in these regions, it does reduce its value by at least a factor of 10 4 .This is sufficient both for computing the energy release due to phase separation and for computing the effect of a neutrally buoyant region on g-mode oscillations.We note that we have also implemented the standard mixing prescription for phase separation which produces uniform profiles in the mixed fluid region.
There is a final simplification we make.Crystallization and phase separation occur in the highly degenerate interiors of WDs, and this degeneracy leads to a decoupling of the thermal and hydrostatic structure.We exploit this as follows.When the model first reaches the temperature at which it begins to crystallize, we tabulate in advance the chemical profile for an arbitrary crystallized mass fraction ranging from 0% to ∼100% of the carbon/oxygen core, taking into account mixing of the fluid layers according to either the standard prescription or our proposed new prescription.We then compute the gravitational/internal energy released by this process, again for an arbitrary value of the crystallized mass fraction.While we do not take into account the effect of the subsequent thermal evolution of the model on the energy release, the error this introduces is quite small.In Montgomery et al. (1999), we showed that this resulted in an error of less than 0.5% in the calculation of the total energy released during phase separation.
With the composition profile and energy release due to phase separation now tabulated as a function of crystallized mass fraction, we simply interpolate at each subsequent time step to find the appropriate values to use for the evolving WD model.

B. THERMOHALINE MIXING
In this section we suggest a possible addition to the Brown et al. (2013) formalism.As a starting point, we first present their approach in section B.1 below, followed by our modification in section B.2. B.1.The Formalism of Brown et al.
The Brown et al. (2013) prescription for thermohaline mixing represents the current state-of-the-art, particularly due to its ability to match the results of numerical simulations over a wide range of parameter values, although these values still do not reach the regimes of most astrophysical objects.These simulations show an initial instability evolving into a saturated, turbulent regime.Since these simulations cannot currently be run for the parameter values found in actual stars, we need a theory or fitting formula that allows us to extrapolate to the cases of interest.
Using the Boussinesq approximation (Spiegel & Veronis 1960), Brown et al. (2013) and Garaud et al. (2015) de-dimensionalize the fluid equations, with length given in units of d = κ T ν kin /N 2 T 1/4 and time in units of d 2 /κ T = Pr 1/2 /N T .This yields the following set of equations: Here, P, T , and µ are the perturbations of the pressure, temperature, and composition on top of an assumed linearly-varying background state, ⃗ u and w are the total fluid velocity and its vertical component, respectively, Pr ≡ ν kin /κ T , the Prandtl number, is the ratio of the kinetic viscosity to the radiative diffusivity, and τ ≡ κ C /κ T is the ratio of diffusivities of the composition to the temperature.R 0 controls the "strength" of the instability, and is defined as where we have used terms defined in equations 5 and 6.Next, Brown et al. (2013) search for linearly growing solutions of the form where the quantities with 0 subscripts are constants, ⃗ l is wavenumber in the x-y plane, and λ is the growth rate of the instability.This leads to an eigenvalue problem with λ given by the solution of a cubic equation (see also Baines & Gill 1969): To find the maximum growth rate, λ m , and its associated wavenumber, l m they derive a new equation by taking a derivative of equation B7 with respect to l and setting dλ/dl = 0.They then simultaneously solve these two equations for λ m and l m .Brown et al. (2013) report their results as Nusselt numbers, which are defined as the total flux divided by microscopic diffusive flux: Next, to find the amplitudes of w and µ, they set the growth rate of a secondary shearing instability equal to that of the primary instability, arriving at the following expression: Thus, Nu µ is a function of Pr, τ , and R 0 .We note that assuming R 0 ∼ O(1) and solving the above equations in the limit τ ≪ Pr ≪ 1 yields In Figure 7, we plot D thrm ∝ Nu µ − 1 as a function of R 0 for the values of τ and Pr given in Table 1.While close to R −1/2 0 , we find the actual scaling to be D thrm ∝ R −0.62 0 .Since the flux F µ is proportional to Nu µ times the µ gradient, which in turn is proportional to 1/R 0 , we find that for our case the flux should scale as To aid in comparing and plotting the different simulations they remap the "fingering regime" of R 0 , 1 < R 0 < τ −1 , onto the unit interval, defining From fits to their numerical simulations, they find C ≈ 7 and are able to provide a remarkably good fit over a wide range of parameters (see Figure 6a).The predicted horizontal size of their "bubbles" is 2π/l m ≈ 10 d, which also matches the simulations well.On the other hand, for the lowest Pr and τ values, there is about a factor of two mis-match between the theoretical and numerically simulated value of Nu µ , and these values are the closest to the regimes of astrophysical interest.

B.2. A New MLT Derivation
In addition to equations B7 -B10 for the eigenvalue λ, we also need relations between the amplitudes of the various perturbations on the RHS of equation B6.In particular, from equation 26 of Brown et al. (2013) we have To estimate the velocity, w 0 , we take the growth rate, λ, times the mixing length, 2π/l, which yields ) where C 2 = (2π) 2 .We note that the only difference between this and the Brown et al. (2013) formula is our prediction for the value of C 2 .We now wish to apply a correction that takes into account turbulent mixing.In a turbulent medium, a fluid element will be continuously mixed with its environment, which will tend to reduce the vertical flux.We can model this as a diffusive process, with an additional diffusion term.For example, we can make the following replacement on the RHS of equation B4: Since D is due to turbulent motions, we estimate it as a velocity times a length scale, i.e., D ∼ ψ λ/l 2 , where ψ is a dimensionless constant we will later adjust.
the addition of a such a term in equation B4 modifies equation B16 so that it becomes Thus, the only modification is that λ → λ(1 + ψ).

(B23)
Adding a term to the right-hand sides of the other equations also results in λ being replaced by λ(1 + ψ).If we define λ 0 ≡ λ(1 + ψ) as the solution of the unmodified eigenvalue equation B7, then we find ) We find that setting ψ = 0.35 in equation B25 results in a good match to the low-Pr, low-τ results.On the other hand, the match to many of the other results is noticeably worse, and the overall fit to the data is not as good as that of Brown et al. (2013).We posit that the above correction should only be applied when the fluid flow is actually turbulent, which should be when the Reynolds number, Re, exceeds some critical value.We use the following ad hoc prescription to calculate ψ: We provide no justification for the Pr dependence of the threshold value η other than it better matches the data.
Figure 6a shows the original fits of Brown et al. (2013) to their simulation data while Figure 6b shows our new fits to these same data.For our fits we also take C = 7, and we choose the value of ψ = 0.35 to match the low-Pr, low-τ results.The threshold η = 2000 Pr was chosen to improve the fits for other cases.For instance, this resulted in an improved fit to the Pr = 1/30, τ = 1/10 case, which is poorly fit by Brown et al. (2013).Of course, our fits should improve since we are increasing the number of fit parameters from 1 to 3. The region R0 > 1 represents the standard thermohaline mixing case, while R0 < 1 is the RT convective region.We see that our new prescription has a mixing efficiency about half as large as that of Brown et al. (2013), for the model WD considered here.
We note a recent paper by Fuentes et al. (2023) presents a mixing length theory for chemical transport that treats in a unified way both the thermohaline and convective regimes.A significant difference with our present work is that they take the mixing length to be of order a pressure scale height, whereas we use a much smaller value of ∼ 10d ≈ 100 cm.While they have made a careful treatment of heat exchange and other diffusive processes, it is not yet clear how their work relates to formal treatments of thermohaline mixing, such as that of Brown et al. (2013).B.3.Spanning the R 0 > 1 and R 0 < 1 Regimes If we consider λ m and Nu µ as a function R 0 (instead of r, see equation B15), then we find that our MLT approach can naturally be extended to the Ledoux convective (RT) regime.In Figure 7 we plot the thermohaline mixing coefficient (in cgs units) as a function of R 0 .We have chosen representative values for Pr, κ T , and κ µ as given in Table 1.As advertised, our prescription is about half as efficient as that of Brown et al. (2013).Each dashed black curve shows a R −0.62 0 scaling of the value of D thrm at R 0 = 1, which is different from but similar to the expected asymptotic scaling of R −1/2 0 (equation B13).Applying this prescription in the RT regime means using L = 2π/l m as the mixing length.While L does increase as R 0 decreases, it is still very much less than the traditional choice of a pressure scale height.It is quite possible, or even likely, that additional instabilities occur in this regime which lead the actual mixing length to be of order a pressure scale height.This would lead to an increase in mixing of several orders of magnitude.Even so, the values shown in Figure 7 should provide a lower bound for the mixing in this regime.A numerical exploration of the transition from thermohaline to RT convection would help clarify this transition.
C. PARAMETER VALUES FOR THE THERMOHALINE CASE Models of thermohaline mixing depend on the microphysical properties of the fluid in question.These can be transport properties (e.g., viscosity and compositional and radiative diffusivities) or thermodynamic properties (e.g., χ T , χ ρ ), and can involve the dimensionless ratios of these quantities (e.g., Pr, τ ).
In Table 1, we list several quantities that are important for the discussion in this paper.We also list approximate values in cgs units for these quantities as obtained from our 1 M ⊙ WD model.In the degenerate interior of a WD, the mean free path of the electrons becomes longer, and this leads to a much higher value for the viscosity than is found in non-degenerate stars.We estimate a Reynolds number for the thermohaline mixing that is ∼ 500 and a Prandtl number, which is the ratio of the kinetic viscosity and the thermal diffusivity, of order Pr ∼ 10 −2 .

Figure 3 .
Figure3.The derived Reynolds number of the mixing for the case shown in Figure2b.We note that our calculation of Re is independent of the mixing length L.

Figure 5 .
Figure5.The flux of C due to different prescriptions for thermohaline mixing compared with the flux required to maintain a chemical profile that is neutrally buoyant (blue curve).We see that theKippenhahn et al. (1980) andBrown et al. (2013) prescriptions can easily accommodate the carbon flux due to phase separation.Only theTraxler et al. (2011) prescription would be too inefficient to produce the required fluxes, presumably leading to the formation of neutrally buoyant profiles.
The authors thank Dean Townsley and Evan Bauer for useful discussions.M.H.M. and B.H.D. acknowledge support from the Wootton Center for Astrophysical Plasma Properties, a U.S. Department of Energy NNSA Stewardship Science Academic Alliance Center of Excellence supported under award numbers DE-NA0003843 and DE-NA0004149, from the United States Department of Energy under grant DE-SC0010623, and from the National Science Foundation under grant No. AST 1707419.M.H.M. acknowledges support from the NASA ADAP program under grant 80NSSC20K0455.

Figure 7 .
Figure7.The thermohaline mixing coefficient (in cgs units) as a function of R0.The region R0 > 1 represents the standard thermohaline mixing case, while R0 < 1 is the RT convective region.We see that our new prescription has a mixing efficiency about half as large as that ofBrown et al. (2013), for the model WD considered here.
The fractional mass of the neutrally mixed plus crystallized region (blue curve) as a function of the crystallized mass fraction.The arrow indicates how a model that is 20% crystallized would have a crystallized plus neutrally buoyant region making it appear 49% crystallized in models that assume uniform mixing.

Table 1 .
Table of Symbol Definitions and Typical Values