A Short Intense Dynamo at the Onset of Crystallization in White Dwarfs

The origin of large magnetic fields (≳106 G) in isolated white dwarfs is not clear. One possible explanation is that crystallization of the star’s core drives compositional convection, which when combined with the star’s rotation, can drive a dynamo. However, whether convection is efficient enough to explain the large intensity of the observed magnetic fields is still under debate. Recent work has shown that convection in cooling white dwarfs spans two regimes: efficient convection at the onset of crystallization, and thermohaline convection during most of the star’s cooling history. Here, we calculate the properties of crystallization-driven convection for cooling models of several white dwarfs of different masses. We combine mixing-length theory with scalings from magnetorotational convection to estimate the typical magnitude of the convective velocity and induced magnetic field for both scenarios. In the thermohaline regime, we find velocities ∼10−6–10−5 cm s−1, with fields restricted to ≲ 100 G. However, when convection is efficient, the flow velocity can reach magnitudes of ∼102–103 cm s−1, with fields of ∼106–108 G, independent of the star’s rotation rate. Thus, dynamos driven at the onset of crystallization could explain the large intensity magnetic fields measured for single white dwarfs.


Introduction
One of the most striking phenomena in astrophysics is the ubiquity of magnetic fields.Like the Earth, stars and stellar remnants such as white dwarfs have a magnetic field.It is known that magnetic fields in a fraction of white dwarfs (WDs) can be more than a million times stronger than that of the Earth ( obs ∼ 10 4 -10 9 G).However, the origin of such large magnetic fields has been a mystery (see Ferrario et al. 2015Ferrario et al. , 2020, for a review), for a review).
Recently, it has been suggested that white dwarfs create magnetic fields through a process similar to the dynamo in the Earth's core.This process involves the crystallization of the white dwarf's interior as it cools down (e.g., Kirzhnits 1960;Abrikosov 1961;Salpeter 1961;van Horn 1968), which causes a separation of its chemical components (Stevenson 1980).The oxygen within the white dwarf preferentially goes into the solid phase, with a corresponding amount of carbon remaining in the liquid phase.This separation results in compositional convection above the crystallized core due to the carbon's buoyancy (Schatzman 1982;Mochkovitch 1983;Isern et al. 1997).With rotation added to the mix, convection can lead to a magnetic dynamo being sustained (Isern et al. 2017;Belloni et al. 2021;Schreiber et al. 2021Schreiber et al. , 2022;;Ginzburg et al. 2022).
Typically, the dynamo's strength is estimated by balancing the energy stored in the magnetic field with the kinetic energy carried by the convective motions  2 /4 ∼  2  (this is known as the equipartition or saturated dynamo hypothesis), where , , and   are the magnitude of the magnetic field, mass density, and convective velocity, respectively (Christensen et al. 2009;Christensen 2010).To estimate   during crystallization, previous studies (e.g., Isern et al. 2017;Ginzburg et al. 2022) have assumed that the kinetic energy flux carried by the convective motions,   , is of the same order as the flux of gravitational energy released from chemical separation,  grav (also ∼ convective heat flux,   ).This leads to estimates for the convective velocities of   ∼ 10 2 -10 6 cm s −1 , which are enough to generate magnetic fields of  ∼ 10 4 -10 9 G, consistent with observations.While some of the assumptions above may be correct, the large thermal conductivity of degenerate electrons in white dwarf interiors can reduce the efficiency of convection (i.e., convection occurs at small Péclet number, Pe, the ratio of thermal diffusion time to convective turnover time).In fact, previous estimations including thermal diffusion gave convective velocities of only 10 −6 -10 −1 cm s −1 , depending on the star's rotation rate (Mochkovitch 1983).
Recently, Fuentes et al. (2023) used mixing-length theory and CIA balance (a balance between the Coriolis, inertial, and buoyancy force; Aurnou et al. 2020) to investigate in detail the transport properties of compositionally-driven convection in white dwarfs.They found that thermal diffusion significantly reduces the convective speed to < 10 −2 cm s −1 , in agreement with the early calculations by Mochkovitch (1983).Further, Fuentes et al. (2023) found that  grav ∼   ≫   , meaning that even though there is enough energy coming out from the star to explain the observed magnetic fields, only a small fraction is in the form of kinetic energy.When combined with the saturated dynamo hypothesis, those much smaller velocities result in magnetic fields of just a few G, in disagreement with observations.Those results have been recently confirmed by Montgomery & Dunlap (2023), who used the transport prescriptions for thermohaline mixing in Brown et al. (2013) to study fluid mixing during phase separation in crystallizing white dwarfs.
The wide range of estimated velocities during crystallizationdriven convection calls for a consistent theory to calculate the flow velocity and the resulting magnetic fields.We identify two ways to improve previous work.First, crystallizationdriven convection spans two regimes during the evolution of cooling white dwarfs: A short-lived regime at the beginning of the crystallization, where the composition flux is large enough so that convection is efficient and Pe ≫ 1, and a second regime of inefficient convection (thermohaline convection) that lasts for most of the cooling history of the star (Castro-Tapia et al. 2024).Second, the large intensity of observed magnetic fields suggest that CIA balance may not be at work in white dwarf interiors, and instead, MAC balance (a balance between the magnetic, buoyancy, and Coriolis force; Stevenson 1979) controls the convective dynamics.
In this paper, we assess the plausibility of crystallizationdriven dynamos by investigating how the MAC balance changes previous estimations for the convective velocities and the resulting magnetic field.In Section 2, we discuss the cooling evolution for a canonical high-mass white dwarf.In Section 3, we derive an analytic expression for the magnetic field induced by crystallization, as a function of different stellar parameters (including Pe).Then, in Section 4 we estimate Pe for both regimes, small Pe or thermohaline convection (relevant for most of the white dwarf's evolution), and the regime of efficient convection or large Pe (important during the early phase of the crystallization).In Section 5 we show our main results, suggesting that dynamos have to be generated at the onset of crystallization.Finally, in Section 6 we summarize and discuss our results.

Cooling history of massive white dwarfs
In this section, we discuss the cooling evolution and crystallization of a 0.9 ⊙ white dwarf.We present a model of a slightly massive white dwarf because strong magnetic fields are mostly observed for WDs of masses > 0.6 ⊙ (Ferrario et al. 2020).Our model is one of those in Castro-Tapia et al. (2024), calculated with the MESA stellar evolution code 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) using a realistic composition profile from stellar evolution.Crystallization starts when the star reaches a luminosity ∼ 10 −3  ⊙ corresponding to a central temperature ≈ 7 × 10 6 K at a central density of 1.7 × 10 6 g cm −3 .The model follows the routine in Bauer (2023) to account for the phase separation of the carbon-oxygen core, and uses the Skye EOS for dense matter in Jermyn et al. (2021), which self-consistently determines the location of the liquid/solid phase transition and the associated latent heat release.More technical details of the model can be found in Castro-Tapia et al. (2024).Figure 1 (top panel) shows the evolution of the core radius and convection zone radius with time.As the star cools, the mass and radius of the core increase with time, while the outer boundary of the convective region moves outward at a different rate.At the onset of crystallization, the convection zone spans < 10% of the star's radius.The C/O profiles at the beginning of crystallization have a well-mixed region that extends further out, up to 0.45 WD (see bottom panel).However, we find that thermal buoyancy stabilizes the unstable composition gradient through most of this region, leaving only a small inner part that is unstable to the Ledoux criterion (≲ 0.1 WD ) initially.
Subsequently, the convective boundary moves out to ≈ 0.45 WD on a timescale ∼ 1-10 Myr.This early evolution of the convection zone is not shown in the models of Isern et al. (2017) and Blatman & Ginzburg (2024) but the extent of the convection zone is similar to our models at  −  crys > 0.1 Gyr.We suppose that this difference would be related to the time resolution of the data points shown in each case.At the same time, the rate at which the convection zone moves outwards in this phase depends on the amount of mixing at the convective boundary, and is therefore sensitive to the prescription used for convective boundary mixing.As a comparison, we also computed models in which the stabilizing effect of thermal buoyancy was ignored, finding that convection then quickly moves out to 0.45 WD (dotted curve in the top panel of Figure 1).The exact evolution of the convection zone in this early phase will require more careful modelling, likely multidimensional numerical simulations of convection.
At later times ≳ 1 Gyr, the convection zone shows a steplike behavior with time.This occurs because of the initial composition profile of the model.The jump in the carbonoxygen abundances temporarily halts the convection zone until it can become carbon-rich enough to push out further (see bottom panel).As a result, the size of the convective region changes significantly along the cooling history of the star.After ∼ 5 Gyr,  out ≈  core ≈ 0.8 WD and the convection zone disappears.At this point, the core contains almost the entire mass of the star.
The efficiency of crystallization-driven convection is set by the strength of the mass flux of light elements (carbon) due to chemical separation at the solid-liquid interface,   .Since   is controlled by the cooling rate of the core, the convective efficiency changes over the cooling history of white dwarfs.As showed by Fuentes et al. (2023) and more recently by Castro-Tapia et al. (2024), a natural dimensionless parameter that measures the strength of   (and the efficiency) is where ∇ ad is the adiabatic temperature gradient, /  is the thermal diffusion time across a pressure scale height   (where   is the thermal diffusivity), and   =   /  is the timescale at which the light elements are released at the solid-liquid interface.
There are two regimes of convection with a rapid transition between them as  increases.For  < 1, the system undergoes thermohaline (inefficient) convection, while for  > 1, the fluid undergoes efficient overturning convection.Figure 2 shows that  ∼ 1-10 2 at the onset of crystallization, and it remains above unity for ∼ 10 Myr.Then, it quickly drops to  ≪ 1 and convection transitions to the thermohaline regime, dominating most of the star's cooling history (∼ 10 Gyr).
We make an analytic model to confirm that the divergence of  towards large values at the onset of crystallization is real 10 3 10 2 10 1 10 0 10 1 t t crys (Gyr)  2024), there is a rapid transition between  > 1, where convection is efficient, and  < 1, the regime of inefficient thermohaline convection.Note that the small fluctuations are numerical artefacts due to the chemical separation which happens in discrete events in the code.For comparison, we include the analytic prediction for  using Mestel's law for cooling white dwarfs (Equation 3, with   = 7 × 10 6 K, and  = −0.1 Gyr) .and not numerical.Assuming () ≈ (4/3)    3 , where   is the central density, and combining hydrostatic balance with the equation of state for degenerate electrons  =   5/3 , where  is a constant, we obtain the density profile () =   (1 −  2 / 2 WD ) 3/2 .The star's interior begins to crystallize when the plasma coupling parameter Γ ∝  1/3 / (the ratio of the electrostatic to thermal energy) exceeds a critical value Γ crit ≈ 175 (van Horn 1968;Potekhin & Chabrier 2000;Bauer et al. 2020;Jermyn et al. 2021).In particular, assuming that Γ crit remains constant during the crystallization, we can write /  = (/  ) 3 , where   is the initial temperature at which the central region of the star begins to crystallize (the external layers of the white dwarf's have a lower density, and therefore reach Γ crit and crystallize at a later time).Then, we can identify the radius of the solid core as a function of temperature where  core is the core temperature.
Then for an incompressible core  core = (4/3)    3 core , we write  core = 4   2 core  core .Taking the time derivative of Equation (2), and using   =  core Δ/4 2 core (where Δ is the carbon enhancement in the liquid phase relative to the solid), Equation 1 can be written as where  = − WD (  /  ) (Δ/) (   /  ∇ ad ) ≈ −0.1 Gyr and it does not vary significantly during the early stages of crystallization.We can see that in this simple model,  is mostly determined by the cooling history of the star's core.Adopting Mestel's cooling law  core ∝  −2/5 (Mestel 1952), we find a reasonable agreement with the output from the MESA model (see dashed curve in Figure 2).As we previously showed in Fuentes et al. (2023),  can be related to the Péclet number Pe (with Pe ∼ 1 separating the two efficiency regimes).Therefore, we must calculate Pe and the convective velocities for both cases (thermohaline and overturning convection) in order to assess the plausibility of crystallization-driven dynamos.

Convective velocities and stellar dynamo
In this section, we combine the MAC balance and scalings for convection in rotating-magnetized flows to derive analytic expressions for the convective velocity during crystallization, and the strength of the induced magnetic field.
Once the dynamo saturates, convective flows are set by a balance between the magnetic, buoyancy, and Coriolis force (MAC balance) where   ≡ / √︁ 4 is the Alfvén speed, ℓ ⊥ is the length scale perpendicular to rotation, and  / is the density contrast in the convective zone.In terms of the composition gradients where   =  ln / ln |  , , ∇  =  ln / ln | ★ is the light element gradient in the star, ℓ is the mixing length, and the critical composition gradient is where  is a constant that depends on the shape of a fluid element (e.g.,  = 9/2 in Kippenhahn et al. 2012).In the limit of large Pe, ∇ ,crit is the gradient needed to be unstable to overturning convection according to the Ledoux criterion.At small Pe, thermal diffusion reduces the effective thermal stratification, so that a smaller ∇ ,crit is needed to drive convection.This is the regime of thermohaline (fingering) convection.For more details, see Fuentes et al. (2023).Davidson (2013) demonstrated that in flows that have small values of the magnetic Prandtl number, Pm ≡ / (where  is the kinematic viscosity, and  is the magnetic diffusivity), the magnetic energy density is determined by the rate of work done by the buoyancy force.This hypothesis and dimensional analysis lead to Typical conditions for the transport properties in the convective region of carbon-oxygen white dwarfs give  ∼ 3×10 −2 cm 2 s −1 ,  ∼ 6×10 −2 cm 2 s −1 , resulting in Pm ∼ 0.5 (Cumming 2002;Isern et al. 2017), suggesting that these scalings should apply.Then combining Equations ( 4) and ( 7), we obtain Since the dynamical length scale for convection in rapidlyrotating-magnetized flows is ℓ ⊥ , let us use ℓ = ℓ ⊥ , so that Pe =   ℓ ⊥ /  .Rewriting Equation ( 8) as shows that the Alfvén speed is entirely determined by the Péclet number, the rotation rate, and the vertical length scale of the fluid motions Equation ( 10) can be re-written in terms of the magnetic field where  rot = 2/Ω, and we have used typical values for cooling white dwarfs.Although Equation ( 11) gives a general expression to evaluate the magnetic field as a function of different stellar parameters, we must first calculate the Péclet number Pe, which can vary dramatically between the thermohaline regime (Pe,  ≪ 1) and the efficient regime (Pe,  ≫ 1).Further, we show later that in the efficient regime, the magnetic field becomes independent of the rotation period  rot , a fact that is supported by observational data.In what follows, we use mixing-length theory to obtain analytic expressions for Pe for each regime.

Regime of small Pe (thermohaline convection)
In a mixture of two elements (e.g.carbon and oxygen), the total composition can be described by only one variable, here chosen to be , the mass fraction of the lighter component (carbon in the context of core crystallization).From mixinglength theory, the mass flux of light elements is where   = ∇  (ℓ/2  ) is the excess composition carried by a fluid element.Equation ( 12) allows us to write the convective velocity as   ≈ (  /∇  ) (2  /ℓ), giving an expression for Pe in terms of the composition flux   , In the limit of Pe ≪ 1, thermal diffusion reduces the effective thermal stratification, i.e., ∇ ≪ ∇ ad , where ∇ =  ln / ln | ★ is the temperature gradient in the star.As a consequence, the critical composition gradient becomes smaller, as we can see by taking the small Pe limit of Equation (6), The reduction in ∇ ,crit allows for composition transport with a composition gradient very close to the critical value, ∇  ≈ ∇ ,crit (see Section 2.1 in Fuentes et al. 2023).Then, Equation ( 13) becomes where we used  = ( therm /  ) (   /  ∇ ad ) and  = 9/2.Typical convection parameters above the crystallization front for a 0.9  ⊙ white dwarf give  ∼ 10 −2 (see Figure 2), then Pe ≈ 0.3, in good agreement with the estimates of Mochkovitch (1983) and Isern et al. (1997).

Regime of large Pe (efficient convection)
Convection in the large Pe regime is different in that the composition gradient ∇  significantly exceeds ∇ ,crit once  > 1 (Fuentes et al. 2023).To find an expression for ∇  , we can balance the Coriolis and buoyancy forces in Equation (4), and rewrite  / in terms of the gradients using Equation ( 5).This gives Using ℓ = ℓ ⊥ , and conveniently introducing   , ∇ ad , and   , this can be rewritten as where are dimensionless parameters that resemble the Taylor and Rayleigh numbers that arise in thermal convection.
For the model in Figure 2,  ∼ 10 (on average) during the efficient regime.Then, using Ra T ∼ 10 26 , and Ta ∼ 10 25 , Equation ( 23) gives Pe ∼ 10 6 , approximately seven orders of magnitude larger than in the thermohaline regime.

Magnetic dynamo at the onset of crystallization
We now evaluate the magnetic field for a set of cooling models of carbon-oxygen white dwarfs with masses  WD ≈ 0.4-1.3⊙ , calculated with MESA2 as described in Section 2. These models were constructed by rescaling a 0.6  ⊙ WD model to the new mass, keeping the same composition profile as a function of relative mass coordinate.We extract the thermodynamic and transport properties above the crystallization front directly from the simulations (see Table 1).For all our models, the thermohaline regime has , Pe ≪ 1 (e.g. ∼ 10 −2 , Pe ≈ 0.3 for the 0.9 ⊙ model, see Section 4), giving   ≲ 10 −5 cm s −1 and  ≲ 100 G for Pe ≲ 1.This shows that thermohaline convection can not explain observed magnetic fields of  obs ∼ 10 4 -10 9 G.We therefore focus on the regime of efficient convection at the onset of crystallization.
It is useful to combine Equations ( 11) and ( 23) to obtain an expression for the magnetic field in the limit of efficient convection in terms of the local stellar properties, Table 1.Thermodynamic and buoyancy properties above the crystallization front for our models of cooling white dwarfs, in the regime of efficient convection.To optimize the construction of the models' grid, with a resolution of 0.1  ⊙ , we use a mass relaxation control in MESA to adjust the WD mass of the settled model in the wd_cool_0.6Mtest suite.As the fast convection occurs on the onset of the crystallization, we increase the time resolution of the simulations before the solid phase appeared, so that the parameters were more accurate at the age  crys , when the crystallization started.For the 1.3  ⊙ model, the outer boundary condition was changed from atm_table='WD_tau_25' to 'photosphere' to help convergence.Note that this is independent of the star's rotation rate.Figure 3 shows the predicted magnetic field strength as a function of white dwarf mass.We see that our calculation for crystallization-driven dynamos using the MAC balance can naturally explain many of the observed values of magnetic field, with the predicted values lying in the range  ∼ 10 6 -10 8 G.The predicted values correlate strongly with white dwarf mass, whereas the observed magnetic fields do not.
It is interesting to note that when separated into H-rich and He-rich, the hydrogen ones seem to follow qualitatively a similar trend to the model (although with large dispersion).Note however that a comparison between the predicted and observed magnetic fields requires modelling the transport of magnetic field from the core where it is generated to the surface (Blatman & Ginzburg 2024).Finally, for the cooling models investigated in this work, Equations ( 8), (10), and (23), give   ∼ 10 2 -10 3 cm s −1 in the efficient regime, increasing with the mass of the star.Comparing the kinetic and magnetic energy densities, we obtain  2 /8 ∼ 10 10 -10 14 erg cm −3 and 1 2  2  ∼ 10 9 -10 14 erg cm −3 , meaning that equipartition is approximately satisfied in the regime of efficient convection.

Discussion
Our results suggest that the crystallization-driven dynamo operates in a short intense phase following the onset of crystallization, before convection transitions to the thermohaline regime.The long Ohmic diffusion timescales in white dwarfs, particularly once the fluid solidifies (e.g.Cumming 2002) means that we can expect the strong field created in these initial stages of crystallization to survive for the remaining lifetime of the white dwarf.
Unlike previous works that rely on energy arguments (Isern et al. 2017) or the CIA balance (e.g., Ginzburg et al. 2022;  24), where the physical parameters above the crystallization front were extracted from a grid of MESA models.The shaded region corresponds to the uncertainties in the predictions due to the dispersion in the measurements of the dimensionless composition flux  at the early phase of crystallization.We include measurements of observed magnetic fields in isolated white dwarfs taken from Bagnulo & Landstreet (2022).We classified the observational data in 4 different groups (H and He rich, and whether they are crystallizing or not based on the age at which the MESA models begin crystallization).Fuentes et al. 2023), our estimations of the convective speed and magnetic fields take into account the correct force balance (Lorentz, Coriolis, and buoyancy force).Our analytic expression for the magnetic field is independent of the rotation rate of the star (see Equation 24), a fact that is supported by obser-vations (e.g., Kawka 2020;Ferrario et al. 2020).For the set of white dwarf cooling models we have used here, the predicted field strength increases with mass (Figure 3) which is not seen in the data, although more work is needed to investigate a more complete set of white dwarf models and to calculate the evolution of the field through the phase of thermohaline convection.The fact that the convection is driven from a small central solid core during the efficient dynamo phase is likely to result in a large angular scale, ie.dipole-dominated, magnetic field geometry.
As pointed out recently by Blatman & Ginzburg (2024), another important problem for crystallization-driven dynamos in white dwarfs is emergence of the magnetic field at the star's surface.Initially, the field is trapped in the convection zone, deep inside the core.Only later on in the evolution of the star, when the convection zone has expanded outward closer to the surface, is the magnetic diffusion time (e.g., Cumming 2002;Saumon et al. 2022) short enough for the field to break out to the surface and be observed.The extent of the convection zone during the efficient convection phase depends on the initial carbon-oxygen abundance profile which is set during stellar evolution and depends on the star's mass (with the convective mantle of massive white dwarfs being closer to the surface, e.g., Bauer 2023).Therefore, the observed strength of the magnetic field will depend mainly on 1) transport of magnetic field between the convection zone and the surface, and 2) how thermohaline convection driven at later times during the cooling history of the star responds to the large magnetic field induced at the onset of crystallization.For this, simulations of magnetized-thermohaline convection, similar to the ones carried out by Harrington & Garaud (2019) and Fraser et al. (2023) will be particularly useful.

Figure 1 .
Figure 1.Top panel: Temporal evolution of the core mass,  core , core radius,  core , and outer radius of the convection zone,  out , once the white dwarf begins to crystallize.Results are shown for a cooling model of a 0.9 ⊙ white dwarf obtained with the MESA stellar evolution code.For  out we show two curves, corresponding to two different ways of determining the extent of the convection zone (the Ledoux criterion or just the unstable composition gradient ignoring the thermal buoyancy).Bottom panel: Chemical profiles of the C/O abundances at different times since the onset of crystallization for the model using the Ledoux criterion shown in the top panel.The profiles at  −  crys = 0.3 Gyr are shown during the regime of thermohaline convection (see text).
Temporal evolution of the dimensionless composition flux at the solid-liquid interface, .Results are shown for a cooling model of a 0.9 ⊙ white dwarf obtained with the MESA stellar evolution code.As discussed in Fuentes et al. (2023) and Castro-Tapia et al. (

Figure 3 .
Figure 3. Intensity of the magnetic field generated at the onset of crystallization, as a function of the white dwarf's mass.The solid line corresponds to predictions from Equation (24), where the physical parameters above the crystallization front were extracted from a grid of MESA models.The shaded region corresponds to the uncertainties in the predictions due to the dispersion in the measurements of the dimensionless composition flux  at the early phase of crystallization.We include measurements of observed magnetic fields in isolated white dwarfs taken fromBagnulo & Landstreet (2022).We classified the observational data in 4 different groups (H and He rich, and whether they are crystallizing or not based on the age at which the MESA models begin crystallization).