© The Authors 2024

1. Introduction

The origin of strong magnetic fields in white dwarf (WD) stars is a hot topic of modern astrophysics (e.g., Bagnulo & Landstreet 2022). In particular, the debate as to whether the magnetic field can be generated in the core of a WD via crystallization-driven convection and dynamo is ongoing, and whether it can subsequently emerge and be observed at the stellar surface is a subject of discussion (e.g., Isern et al. 2017; Ginzburg et al. 2022; Fuentes et al. 2023; Montgomery & Dunlap 2024; Blatman & Ginzburg 2024, and references therein). Most recently, Fuentes et al. (2024) have proposed a short-lived (≲10 Myr) phase of efficient composition-driven convection characterized by a high upward light element flux and rapid fluid motions that operates at the initial stage of C/O WD crystallization. Lanza et al. (2024) discuss a powerful dynamo driven by the rapid upward flow of ²²Ne-depleted solid flakes in the process of ²²Ne distillation prior to core crystallization in C/O/²²Ne WDs.

The present paper focuses on crystallizing C/O WDs. We studied the thermal evolution of matter in the inner core region at the freezing onset by numerically solving the heat equation on a fine spatial and temporal grid and by including a detailed treatment of the latent heat release. As a byproduct, we obtain an alternative estimate of the light element flux induced by the C/O mixture crystallization.

2. Model

When the freezing of a WD core begins, the very central region reaches its crystallization temperature and then cannot cool further¹ unless all of its latent heat has been released and removed. The release of the latent heat is controlled by the ability of the system to carry it away, that is, by the available thermal gradients and heat fluxes. More outward layers with time also reach their crystallization temperature and start releasing latent heat, that is, the region that has begun crystallization gradually expands. At some point, the central part completely exhausts its latent heat and resumes cooling according to the standard heat equation. This means that, in general, there are three regions when crystallization in a WD begins. In the inner and outer regions, matter cools according to the standard heat equation and is in a crystal and liquid state, respectively. The intermediate region represents a reservoir of heat, where matter is partially solid and partially liquid. Its temperature is equal to the local crystallization temperature². To model the intermediate region, one simply has to keep track of the balance between the total amount of available heat and the incoming and outgoing heat fluxes.

Since we are mostly interested in a qualitative picture of the freezing onset, our model is significantly simplified. Following Fuentes et al. (2024), we assumed a specific form of the mass density, ρ, dependence on the radial coordinate, r,

$$\begin{array}{c}\hfill \rho ={\rho }_{\mathrm{c}}{(1-\frac{r^2}{R_{\mathrm{WD}}^2})}^{3/2},\end{array}$$(1)

where ρ_c = 1.86 × 10⁷ g/cc is the central density, and R_WD = 6371 km is the WD radius. The specific values are taken from the evolutionary model³ of a WD with a thick hydrogen envelope and a mass of M_WD = 0.9 M_⊙ near the crystallization onset (Bédard et al. 2020, the reference model).

We further assumed equal proportions of C and O by mass but neglected any effects of phase separation and element redistribution, including the gravitational energy release, on the thermal evolution. In reality, the liquid around newly forming crystals becomes depleted of the heavier element (oxygen) and enriched with the lighter element (carbon). This creates composition-driven instability conditions that were discussed by, for example, Fuentes et al. (2023) and Montgomery & Dunlap (2024). They cause the upward light element flux, F_X, and an opposite heavy element flux, which replenishes the oxygen content of the liquid. Since we limited our focus to the inner core region and to the initial stages of crystallization, we do not describe these effects in detail, except for estimating the flux F_X and assuming a fixed liquid composition. We expect that including these effects, particularly the additional heat that is produced by mixing, which also has to be removed, would only strengthen our conclusions. Furthermore, we neglected any effects associated with the distillation of ²²Ne (Isern et al. 1991; Blouin et al. 2021) or any other heavier species.

The crystallization (melting) temperature, T_m, was calculated as

$$\begin{array}{c}\hfill T_{\mathrm{m}}=\frac{Z_1^{5/3}{e}^2}{a_{\mathrm{e}}{k}_{\mathrm{B}}{\mathrm{\Gamma }}_{1\mathrm{m}}}{\xi }_{\mathrm{CO}}(3/7),\end{array}$$(2)

where Z₁ = 6 is the carbon charge number, a_e = (4πn_e/3)^−1/3, n_e is the electron density, Γ_1m = 175.6 is the Coulomb coupling parameter of a classic one-component plasma at crystallization (Baiko & Chugunov 2022), and ξ_CO(3/7)=1.05 is a correction at the value of the oxygen number fraction, x_O = 3/7, based on the phase diagram of the C/O mixture (Blouin & Daligault 2021). Furthermore, we assumed: a constant heat capacity of C = 2.7k_B per ion, which is a reasonable approximation for mildly quantum liquids and solids (in our system, the ion plasma temperature is ∼3 times greater⁴ than the actual temperature) in the vicinity of the phase transition (e.g., Baiko & Chugunov 2022); a constant thermal conductivity of κ = 1.2 × 10¹⁶ erg cm⁻¹ s⁻¹ K⁻¹ based on the electron-ion collision frequency calculation, fit, and extrapolation to mixtures (Potekhin et al. 1999); and the latent heat of crystallization for the classic equimass C/O mixture, λ = 0.67 k_BT_m per ion (Baiko 2023).

A realistic study of the crystallization onset requires a description of the local latent heat release and removal, which in turn requires a non-isothermal approach. In our model, the initial temperature at the stellar center is taken to be 0.03% greater than the melting temperature, T_m(0). The initial temperature profile, T_ini(r), is based on the following scaling:

$$\begin{array}{cc}\hfill T_{\mathrm{ini}}(r)& =1.0003T_{\mathrm{m}}(0)-0.5K{r}^2,\hfill \\ \hfill K& =\frac{\rho C}{3\kappa ⟨A⟩m_{\mathrm{u}}}\left|\frac{\mathrm{d}T}{\mathrm{d}t}\right|,\hfill \end{array}$$(3)

where ⟨A⟩ is the average mass number of the mixture, m_u is the atomic mass unit, and dT/dt is a time derivative of the temperature.

We analyzed the thermal evolution of matter only out to r = r₀ with r₀ = 3000 km to avoid the need to describe the complicated physics of the outer, partially ionized and partially degenerate layers. Consequently, a boundary condition at r₀ is required, and for this, we specified the time derivative of the temperature, dT(r₀)/dt = α₀(t). In principle, the true form of α₀(t) can only be obtained from a full-star evolutionary modeling with a high-level code that accounts for the effects investigated in the present paper. In the absence of this, we can first assert that T(r₀) decreases monotonously with t, and we can then take advantage of the fact that we are interested in a relatively brief time interval, Δt <  100 Myr since the crystallization onset, which covers the purported efficient convection stage. For this reason, the time derivative was assumed to be constant, α₀ = −1.5 × 10⁻¹⁰ K/s. This value of α₀ is close to the average time derivative of the central temperature in the reference model shortly before the crystallization onset. On the other hand, in our model this value of α₀ results in an average time derivative of the central temperature after the crystallization onset of −0.83 × 10⁻¹⁰ K/s, which is practically the same as the value in the reference model. In addition, we tried several other values of α₀ that are within a factor of ∼2 of the above fiducial value. This yielded a qualitatively similar picture of the temperature variation. We therefore do not expect that implementing a weak decrease in |α₀(t)| over the considered Δt would affect our main conclusions. We also set dT/dt to α₀ in Eq. (3).

3. Results

The result of the modeling is shown in Fig. 1. We display the temperature profile as a function of r, and the different curves correspond to different time stamps. The dotted red curve demonstrates the initial profile. The four solid green curves between it and the long-dashed blue profile show the evolution before the latent heat in the center is exhausted: the upper solid curve is separated by 2 Myr from the initial curve, and the other three curves add 8 Myr of evolution each. On the left side, these solid curves approach the thick short-dashed orange line, which is the profile of the crystallization temperature, and then smoothly merge onto it. At r = 0, they are equal to T_m(0). In this case, there are therefore only two regions: the inner region is at the local crystallization temperature, and the outer region represents the cooling liquid. In the inset, the top left corner of the graph is shown in more detail. The five solid green profiles to the right of the crystallization temperature are separated from each other and from the initial dotted red profile by 400 kyr.

Fig. 1. Internal temperature profile of the cooling WD model at various times. We show the initial profile (dotted red), the profile shortly after the formation of the fully crystallized core (long-dashed blue), three profiles at 45.6 Myr since the beginning of the modeling (dot-dashed green for the default latent heat value, dot-dashed brown for twice lower and zero latent heat values), and profiles at other times (solid green; see the main text for details). The crystallization temperature is shown with thick short orange dashes.

The long-dashed blue curve corresponds to 27.6 Myr since the beginning of the simulation, slightly later than the moment at which the crystallization of the innermost core region is complete, that is, when all its latent heat is removed. Consequently, the solid green curves below it contain three segments that correspond to the three regions described in the first paragraph of Sect. 2. The segments, describing the cooling crystal with progressively lower central temperatures, are seen to the left of the crystallization profile. The mutual time separation between these curves is 8 Myr (the dot-dashed green curve is also part of this family; see below), and the upper curve is at 2 Myr after the long-dashed profile. This regime is characterized by an expansion of the region that has reached T_m(r) and of the region that has exhausted its latent heat and cooled to below T_m(r). The radial extent of the intermediate region shrinks, and by the last curve, which displays the temperature profile at 69.6 Myr, its size becomes smaller than the spatial grid. The five solid green curves to the left of the crystallization temperature profile in the inset are separated from each other and from the long-dashed blue curve by 400 kyr.

For illustration, we show with dot-dashed curves the temperature profile at 45.6 Myr for three sequences with different latent heat values (assuming the same initial profile and boundary condition): the upper green curve corresponds to the standard λ value, and the brown curves in the middle and below demonstrate the cases with λ = 0.33 k_BT_m and zero latent heat, respectively. For the latter, there should be no difference between the liquid and solid phases in our code. A smooth curve that rapidly cools in the center is in fact obtained.

In Fig. 2 we show the radii of the regions that have started and completed crystallization, denoted as r_m and r_c, respectively, with dot-dashed red and solid blue curves as functions of time. As in Fig. 1, the difference between the initial time stamps for the two processes is ∼27 Myr. At the beginning, both dependences have an infinite time derivative that qualitatively agrees with the analytic model of Fuentes et al. (2024). Eventually, the two curves merge (to within the spatial grid size).

Fig. 2. Radius of the region that has started crystallization (dot-dashed red) and of the fully crystallized core (solid blue) vs. time since the beginning of the modeling.

We next estimated the flux of the light elements, F_X, that results from the phase separation that accompanies crystallization (in our case, this is the flux of carbon ions). Fuentes et al. (2024) link the flux with the velocity of the crystallization front. Specifically, if the front moves from r − dr to r over a time interval dt, then the entire spherical shell of thickness dr goes from the liquid to the solid phase in their model. Based on this, the mass of released light elements that then cross a sphere of radius r over the period dt (by convective motions) can be calculated. This quantity was found to be divergent at the initial stages of crystallization.

However, as our solution shows, this picture does not apply at the crystallization onset. The light elements are released rather slowly within the spherical shell with an outer radius r_m and an inner radius r_c (r_c = 0 for the first ∼27 Myr). The mass of released light elements over an interval dt is proportional to the total mass that has solidified. This in turn is proportional to the amount of latent heat that is carried away from this region. When the convection is sufficiently rapid, then the light elements can cross the sphere at r_m over the period dt. This provides an upper estimate of the convective flux.

We calculated the flux of light elements using both prescriptions and show the results in Fig. 3. Based on the phase diagram of Blouin & Daligault (2021), we set the difference in the carbon mass fraction between the liquid and the solid to ΔX = 0.11. The dot-dashed blue curve assumes that the flux is determined by the fully crystallized core growth rate. In agreement with the analysis of Fuentes et al. (2024), this curve is divergent. The solid green curve⁵ takes the gradual release of light elements into account, which is controlled by the latent heat removal that occurred over a period of ∼27 Myr before the start of the fully crystallized core formation. As soon as the central region reaches the crystallization temperature, the flux rapidly grows, starting from zero. Then it is finite everywhere, grows slowly, and is more than twice lower than the previous estimate, even 10 Myr after the appearance of the fully solidified core. At the end of our modeling, r_c catches up with r_m, and the two calculations of the light element flux coincide, as expected.

Fig. 3. Light element flux vs. time. The dot-dashed blue and solid green curves show fluxes based on the fully crystallized core growth rate and on the gradual release of light elements controlled by the latent heat removal, respectively.

4. Conclusion

We have modeled the thermal evolution of the internal region of a 0.9 M_⊙ C/O WD over a period of ∼70 Myr since the beginning of ion plasma crystallization in its center. We observed the formation of two spherical shells. The outer shell surrounds the region in which crystallization has begun, and the inner shell bounds the fully solidified core, which has exhausted its latent heat. The region between the shells is partially liquid and partially solid. It gradually emits the latent heat generated by the crystallization, which is carried away along the existing thermal gradients, and it also releases light elements in the process of element redistribution that accompanies the mixture solidification. The inner shell first appears ∼27 Myr after the appearance of the outer shell, whose radius reaches ∼1000 km at this moment. At ∼70 Myr, the inner shell catches up with the outer shell, and the crystallization process adopts the familiar motion of a 2D front.

Assuming that all released light elements cross the outer shell, we estimated their flux induced by the mixture crystallization. Until the time when the two shells merge, the result we obtain is significantly lower than the estimate based on the model of Fuentes et al. (2024). In particular, according to Fig. 2 of that work, the efficient convection driven by the light element flux is expected to cease after ∼3 Myr since the formation of the fully solidified core. However, at this moment, the actual flux of the light elements is about five times lower than the Fuentes et al. (2024) prediction and is even lower at earlier times (see Fig. 3). This means that the criterion for efficient convection is not fulfilled at the early stages of WD crystallization, and consequently, such a phenomenon probably does not occur.

Acknowledgments

The author is deeply grateful to D. G. Yakovlev for discussions and to expert referee, Mike Montgomery, for valuable comments. This work has been supported by the Russian Science Foundation grant 24-12-00320.

References

All Figures

Fig. 1. Internal temperature profile of the cooling WD model at various times. We show the initial profile (dotted red), the profile shortly after the formation of the fully crystallized core (long-dashed blue), three profiles at 45.6 Myr since the beginning of the modeling (dot-dashed green for the default latent heat value, dot-dashed brown for twice lower and zero latent heat values), and profiles at other times (solid green; see the main text for details). The crystallization temperature is shown with thick short orange dashes.

In the text

	Fig. 2. Radius of the region that has started crystallization (dot-dashed red) and of the fully crystallized core (solid blue) vs. time since the beginning of the modeling.
In the text

	Fig. 3. Light element flux vs. time. The dot-dashed blue and solid green curves show fluxes based on the fully crystallized core growth rate and on the gradual release of light elements controlled by the latent heat removal, respectively.
In the text