²²Ne Phase Separation as a Solution to the Ultramassive White Dwarf Cooling Anomaly

White dwarf evolution is often depicted as a simple, uneventful cooling process, implying that they naturally lead to accurate age determinations. White dwarfs are indeed very useful cosmic clocks. They have been used to determine the ages of individual stellar populations (Hansen et al. 2007) and to reconstruct the star formation history of the Milky Way (Tremblay et al. 2014; Kilic et al. 2017; Fantin et al. 2019). When paired with a stellar companion, they also become useful benchmarks for harder-to-model objects (e.g., M, L, and T dwarfs, Lam et al. 2020; Meisner et al. 2020). They could even prove useful to track the galactic evolution of lithium (Kaiser et al. 2021).

However, some aspects of white dwarf cooling remain poorly understood, meaning that the theoretical cooling tracks routinely used to infer white dwarf ages may not be as accurate as they ought to be. The existence of those inaccuracies has been clearly demonstrated during the last few years thanks to the Gaia Data Release 2 (Gaia Collaboration et al. 2016, 2018a). Cheng et al. (2019) revealed the most salient example of the shortcomings of current evolution models by showing that some ultramassive white dwarfs (M_⋆ ≳ 1.05 M_⊙) experience an additional cooling delay of ∼8 Gyr compared to theoretical predictions. Those objects are found in a region of the Gaia color–magnitude diagram (CMD) known as the Q branch (Gaia Collaboration et al. 2018b), which corresponds to an overdensity of objects that coincides with the high-mass tail of core crystallization (Tremblay et al. 2019).

White dwarfs experiencing this additional cooling delay most likely have C/O cores, as the C/O crystallization sequence (and not the O/Ne crystallization sequence) is consistent with the location of the observed overdensity in the Gaia CMD (Bauer et al. 2020, Figure 3). In apparent contradiction with this observation, standard single-star evolution models predict the formation of O/Ne cores in ultramassive white dwarfs (Siess 2007), and ultramassive white dwarfs formed from the merger of two C/O white dwarfs are also predicted to have O/Ne cores (Schwab 2021). However, Althaus et al. (2021) recently proposed two single-star evolution scenarios that can lead to the formation of ultramassive C/O white dwarfs.

Assuming a mass fraction of ²²Ne X(²²Ne) = 0.02, ¹ Cheng et al. (2019) note that, if released during crystallization, the gravitational energy of ²²Ne stored in ultramassive C/O white dwarfs is sufficient to generate a ∼6–9 Gyr cooling delay. Such a large effect is possible owing to the neutron-rich nature of ²²Ne (A > 2Z). The most commonly discussed mechanism to release this potential energy is the gravitational settling of ²²Ne in the liquid phase (Bildsten & Hall 2001). However, current models predict that simple gravitational settling is not efficient enough to give rise to multi-Gyr cooling delays (e.g., Bauer et al. 2020), in part because crystallization of the C/O core strongly inhibits the extent of the diffusion of ²²Ne in the liquid phase.

Camisassa et al. (2020) proposed that the objects undergoing this extra delay are extremely rich in ²²Ne, with a ²²Ne mass fraction X(²²Ne) = 0.06 instead of the X(²²Ne) = 0.014 value expected for white dwarfs issued from the single-star evolution of solar-metallicity progenitors (Cheng et al. 2019). While quadrupling X(²²Ne) can indeed produce a multi-gigayear cooling delay, more work is needed to explain how ultramassive white dwarfs could acquire such a high ²²Ne abundance. Simulations by Staff et al. (2012) show that mergers of He and C/O white dwarfs can lead to a sizeable ²²Ne enrichment. Yet, the stars produced in those merger simulations had a mass of only 0.9 M_⊙, meaning that those results may not be directly applicable to ultramassive white dwarfs. Another avenue, explored by Bauer et al. (2020), is to increase the rate of ²²Ne gravitational settling through the formation of solid clusters containing a few thousand ²²Ne ions. However, subsequent molecular dynamics simulations have shown that the formation of such clusters cannot take place at the low X(²²Ne) found in C/O white dwarfs (Caplan et al. 2020). Other avenues unrelated to ²²Ne sedimentation have also been explored, but cannot explain the missing multi-gigayear delay (Horowitz 2020).

In this Letter, we accurately calculate the melting curve of C/O/Ne mixtures to show that the gravitational energy of ²²Ne stored in ultramassive white dwarfs can be efficiently released through a phase separation process. Assuming a modest ²²Ne enrichment (half of that assumed by Camisassa et al. 2020), this mechanism leads to multi-gigayear cooling delays that can largely explain the ultramassive white dwarf cooling anomaly. We also argue that this process can solve a separate, smaller cooling delay problem that affects standard-composition white dwarfs.

The possibility of ²²Ne phase separation ² and its effects on the cooling of C/O white dwarfs were first discussed by Isern et al. (1991). The general idea is that depending on the composition of the C/O/Ne plasma and the exact shape of the C/O/Ne phase diagram, the solid crystals that are formed when the crystallization temperature is reached can be depleted in ²²Ne with respect to the liquid mixture. If the ²²Ne depletion is large enough, those crystals are lighter than the surrounding liquid and float upward, away from the crystallization front. The rising crystals eventually melt in lower density regions where their constituent ions are mixed via Rayleigh–Taylor instabilities. This is analogous to a distillation ³ process and it gradually displaces ²²Ne-rich liquid downward. The ²²Ne abundance in the liquid at the liquid–solid interface thereby increases until it reaches a critical value for which the C/O/Ne phase diagram predicts that there is no more phase separation (i.e., when the liquidus and solidus meet). The plasma then freezes at this constant composition, forming solid ²²Ne-rich layers. The crystallization of the remaining liquid C/O mixture (now free of ²²Ne) then continues as usual.

The amount of gravitational energy released in this process critically depends on the exact shape of the C/O/Ne phase diagram. If the distillation process starts immediately at the beginning of the crystallization of the C/O core (i.e., when the core is still fully liquid), then a very large cooling delay is possible. In this scenario, favored by the exploratory calculations of Isern et al. (1991) and Segretain et al. (1994), all the ²²Ne from the liquid core is transported to the center of the white dwarf where it forms a compact ²²Ne-rich central core. It is also possible, as suggested by the C/O/Ne phase diagram of Segretain (1996), that the distillation process only starts later, when a fraction of the core is already solidified. A ²²Ne-rich shell is then formed around a central core with a nearly unperturbed ²²Ne concentration. In this scenario, a smaller amount of ²²Ne is transported downward and not as deep in the gravitational well of the star, with a more modest energy release.

Caplan et al. (2020) recently used the semianalytic approach of Medin & Cumming (2010) to map the C/O/Ne phase diagram. ⁴ While they do not discuss this possibility, we note that their Figure 4 confirms that the solid phase can be depleted in ²²Ne when the ²²Ne concentration remains small in the liquid phase, enabling the distillation process described above. We also remark that their phase diagram indicates that ²²Ne distillation will continue (the ²²Ne concentration in the liquid will remain higher than that in the solid, ${x}_{\mathrm{Ne}}^{{\ell }}\gt {x}_{\mathrm{Ne}}^{s}$) ⁵ until (x_C, x_O, x_Ne) ≃ (0.8, 0.0, 0.2). This composition corresponds to where the liquidus and solidus meet in the two-component C/Ne phase diagram (see also Figure 4 of Medin & Cumming 2010 and Figure 5 of Ogata et al. 1993).

While past calculations of the C/O/Ne phase diagram support the existence of ²²Ne distillation, they are not precise enough to allow an accurate quantitative assessment of the impact of this process on white dwarf cooling. In particular, a high-resolution version of the C/O/Ne phase diagram at small ²²Ne concentrations is required. The exact shape of the phase diagram in this region is critical, as it dictates even the qualitative outcome of the distillation process (the formation of a ²²Ne-rich central core versus a ²²Ne-rich shell closer to the surface).

To obtain a high-resolution, high-precision version of the C/O/Ne phase diagram at small ²²Ne concentrations, we turn to the Clapeyron integration technique that we have recently developed to map the phase diagrams of dense plasmas (Blouin & Daligault 2021) and applied to the two-component C/O plasma (Blouin et al. 2020). Briefly, this method consists of directly integrating the liquid–solid coexistence line using the appropriate Clapeyron equation for phase transitions at constant temperature and pressure (here, Equation (A2) of Blouin & Daligault 2021), which we evaluate using Monte Carlo simulations where the full electron–ion plasma is considered. All our simulations use N = 686 ions, which is enough to mitigate finite-size effects even at very small concentrations (x_Ne ∼ 0.002). We verified this point by performing simulations with up to N = 4000 ions. Each simulation is executed for 7 × 10⁶ Monte Carlo iterations, of which the last 5 × 10⁶ are used to evaluate the average thermodynamics quantities needed to integrate the Clapeyron equation. We assume a fixed P = 10²⁴ dyn cm⁻² pressure, a typical value for white dwarf cores. We integrate the melting line at constant temperature T, meaning that many distinct integrations are needed to fully map the C/O/Ne phase diagram in the three-dimensional composition–temperature space (we use 11 different temperatures). It is useful to express the temperature as

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{\rm{C}}}=\displaystyle \frac{{e}^{2}}{{a}_{e}{k}_{B}T}{Z}_{{\rm{C}}}^{5/3},\end{eqnarray} \tag{ 1 }$$

where e is the elementary charge, ${a}_{e}={\left(3/4\pi {n}_{e}\right)}^{1/3}$ with n_e the electron density, k_B is Boltzmann's constant, and Z_C = 6. At constant P, n_e is essentially constant in dense degenerate plasmas and Γ_C is therefore a useful dimensionless measure of T. For reference, a pure C plasma solidifies at Γ_C ≃ 175. The 11 temperatures we simulate are equivalent to Γ_C ranging from 132 to 182. Those values in turn correspond to the melting temperatures, Γ_C,m , of C/O plasmas with ${x}_{{\rm{O}}}^{{\ell }}$ ranging from 0.28 to 0.79.

Figure 1 shows the result of a low-resolution integration of the phase diagram at Γ_C = 179.8 (the melting temperature of a C/O plasma with ${x}_{{\rm{O}}}^{{\ell }}=0.33$) between ${x}_{\mathrm{Ne}}^{{\ell }}=0.0$ and ${x}_{\mathrm{Ne}}^{{\ell }}=0.2$. As with the two-component C/O phase diagram, we find that the solid is enriched in O compared to the liquid. As for Ne, the solid can either be enriched or depleted with respect to the liquid depending on the Ne abundance in the liquid. For the temperature used in Figure 1, the solid is depleted in Ne as long as ${x}_{\mathrm{Ne}}^{{\ell }}\lesssim 0.15$. It is in this regime that ²²Ne distillation is possible.

Zoom In Zoom Out Reset image size

To precisely determine when the distillation process begins, we performed a set of high-resolution integrations of the phase diagram between ${x}_{\mathrm{Ne}}^{{\ell }}=0.00$ and ${x}_{\mathrm{Ne}}^{{\ell }}=0.05$. For each point along the coexistence line, we evaluate the mass densities ρ of both phases. Figure 2 shows how the solid–liquid density difference varies as a function of ${x}_{\mathrm{Ne}}^{{\ell }}$ for 3 of the 11 temperatures for which we have integrated the C/O/Ne phase diagram. Temperatures and ²²Ne concentrations that are below the ρ^s = ρ^ℓ line correspond to conditions where ²²Ne distillation operates. This regime is reached at low melting temperatures (corresponding to a high Γ_C,m and a low O concentration) and/or high ²²Ne concentrations. Note that even though we found that the solid is always depleted in ²²Ne in the ${x}_{\mathrm{Ne}}^{{\ell }}\to 0$ limit, its density is not always smaller than that of the liquid as its O enrichment must also be considered.

Zoom In Zoom Out Reset image size

Figure 3 translates the results of Figure 2 in the $({{\rm{\Gamma }}}_{{\rm{C}}},{x}_{\mathrm{Ne}}^{{\ell }})$ plane for all 11 simulated temperatures. A key finding is that both scenarios described in Section 2 (the formation of a ²²Ne-rich central core or shell) are possible depending on the white dwarf's initial composition.

Zoom In Zoom Out Reset image size

If the central composition of the white dwarf is rich enough in ²²Ne (scenario (a) in Figure 3), distillation will start directly at the onset of crystallization when the melting temperature Γ_C,m is reached in the central layers. The core will remain entirely liquid until the distillation process described in Section 2 increases the ²²Ne concentration in the central layers to the critical x_Ne = 0.2 value (the solid crystals keep floating upward, so no solid core can be formed). This ²²Ne enrichment process is symbolized by the blue arrow in Figure 3. The central layers will then freeze at the composition (x_C, x_O, x_Ne) = (0.8, 0.0, 0.2), which is reached at Γ_C,m = 208. After that, the crystallization of the remaining C/O core will proceed as usual. This will lead to the formation of a ²²Ne-rich C/Ne solid core surrounded by a C/O mantle completely depleted of ²²Ne. This final state can also be attained if the central mixture is O poor, in which case the melting temperature is lower and distillation can start right at the beginning of the crystallization process even with a small ²²Ne concentration. This corresponds to scenario (b) in Figure 3.

If a more standard composition is assumed (e.g., scenario (c), X(O) = 0.60 and X(²²Ne) = 0.014, or x_O = 0.53 and x_Ne = 0.009), then no ²²Ne distillation can initially occur and crystallization takes place as in the case of a two-component C/O plasma without any significant change to the ²²Ne distribution. As the crystallization front progresses outward, the liquid is gradually depleted in O due to C/O phase separation and it freezes at increasingly lower temperatures (Horowitz et al. 2010; Althaus et al. 2012; Blouin et al. 2020). Eventually, Γ_C,m crosses the boundary that delimits the distillation regime and ²²Ne distillation can start (this corresponds to the tip of the orange arrow in Figure 3). The progression of the crystallization front is halted until the distillation process increases the ²²Ne abundance in the liquid at the liquid–solid interface to x_Ne = 0.2. The ²²Ne-rich layers surrounding the solid core then freeze, and the crystallization of the remaining C/O layers continues following the two-component C/O phase diagram. This leads to a final structure where a C/O/Ne central core is surrounded by a ²²Ne-rich C/Ne shell, which is itself surrounded by a C/O outer shell completely depleted of ²²Ne.

How much does this process affect white dwarf cooling? To answer this question, we generate white dwarf structures using STELUM (Bédard et al. 2020; A. Bédard et al. 2021, in preparation) and compute their binding energies

$$\begin{eqnarray}&&B=-\left({\int }_{0}^{{M}_{\star }}{udm}-{\int }_{0}^{{M}_{\star }}G\displaystyle \frac{m}{r}{dm}\right),\end{eqnarray} \tag{ 2 }$$

where u is the internal energy per unit mass, G is the gravitational constant, and the integrals are performed over the whole star. The net energy released by the distillation process corresponds to the difference in binding energy, ΔB, between the white dwarf structure at the onset of ²²Ne phase separation and the structure obtained once it ends.

For the initial structure, we assume that ²²Ne is homogeneously mixed throughout the C/O core. This is not strictly correct as ²²Ne gravitational settling in the liquid phase will create a ²²Ne concentration gradient, but this effect is small (e.g., Figures 2 and 3 of García-Berro et al. 2008) and can be neglected for our purpose. For simplicity, we also assume that the C/O abundance ratio is constant throughout the core. For the final structure, we place all the ²²Ne in a central core of mass M_⋆ X(²²Ne)/X^a (²²Ne), where M_⋆ is the white dwarf's mass, X(²²Ne) is its overall ²²Ne mass fraction, and X^a (²²Ne) = 0.3 is the ²²Ne mass fraction in the central core (i.e., x_Ne = 0.2). If a ²²Ne-rich shell is formed instead, then its mass is given by ${M}_{\star }^{{\ell }}X{(}^{22}\mathrm{Ne})/{X}^{a}{(}^{22}\mathrm{Ne})$, where ${M}_{\star }^{{\ell }}$ is the mass of the fraction of the core that is still liquid when the distillation process starts. Table 1 gives the resulting ΔB for different scenarios. We also estimate the resulting cooling delay by dividing this ΔB by the average luminosity of the star over the distillation process. ⁶

The first portion of Table 1 gives the impact of ²²Ne phase separation on ultramassive ²²Ne-rich white dwarfs that will form a ²²Ne-rich central core as in scenario (a) of Figure 3. For those models, we use the X(²²Ne) = 0.035 value discussed in Bauer et al. (2020) for white dwarfs that descend from stars formed in α-rich environments. However, our conclusions are not limited to this particular scenario; any other scenario that leads to a moderate ²²Ne enrichment can be envisaged (e.g., the merger of two white dwarfs; Blouin et al. 2020; Camisassa et al. 2020). We find that ²²Ne phase separation induces an additional 4–9 Gyr cooling delay for those objects, which explains most of the ultramassive white dwarf cooling anomaly. An abundance of X(²²Ne) = 0.035 is expected to be rare in white dwarfs, which is consistent with the fact that only ∼6% of ultramassive white dwarfs are affected by the cooling anomaly (Cheng et al. 2019). Depending on the mass and composition of the white dwarf, the delay induced by ²²Ne phase separation can be short of the 8 Gyr delay inferred by Cheng et al. (2019), but there are additional factors to consider:

• 1.  
  The ²²Ne mass fraction for the delayed ultramassive white dwarfs may well be higher than the X(²²Ne) = 0.035 value assumed here.
• 2.  
  As explained in Bauer et al. (2020), the missing cooling delay can be reduced to ∼6 Gyr by changing some assumptions in Cheng et al.'s analysis, such as the thick disk age, the age–velocity dispersion relation, and the star formation history.
• 3.  
  The C/O white dwarf cooling tracks used by Cheng et al. (2019) for M_⋆ < 1.10 M_⊙ do not include C/O phase separation (Fontaine et al. 2001), meaning that they underestimate the time those objects stay on the Q branch by up to ∼1 Gyr (Blouin et al. 2020).
• 4.  
  Cheng et al. (2019) use O/Ne white dwarf cooling tracks (Camisassa et al. 2019) for M_⋆ ≥ 1.10 M_⊙, but it is increasingly clear that the objects that form the delayed population have C/O cores (Bauer et al. 2020; Camisassa et al. 2020). As pointed out by Bauer et al. (2020), O/Ne white dwarfs crystallize much earlier in their evolution than C/O white dwarfs (i.e., at higher luminosities), meaning that the latent heat and gravitational energy released during crystallization induce a shorter cooling delay and that ²²Ne diffusion has less time to operate. From Figure 2 of Camisassa et al. (2020), we infer that this led Cheng et al. (2019) to overestimate by ∼1.5 Gyr the cooling delay experienced by M_⋆ ≥ 1.10 M_⊙ objects assuming they have C/O cores.

Overall, the cooling delay provided by ²²Ne phase separation appears to be the key piece of missing physics required to explain the ultramassive white dwarf cooling anomaly. Note that our results provide further support for the idea that the delayed population identified by Cheng et al. (2019) have C/O cores, as ²²Ne distillation cannot take place in O/Ne cores. Neglecting quantum effects in the ionic degrees of freedom, the phase diagram does not discriminate between ²⁰Ne and ²²Ne. The total ²⁰Ne+²²Ne abundance being higher than that of the azeotrope (Figure 1 of Camisassa et al. 2019), the solid will be enriched in ²⁰Ne and ²²Ne, making it denser than the liquid. A sedimentation process analogous to that occurring in C/O white dwarfs will take place.

Any solution to the ultramassive cooling anomaly must also preserve the existing agreement between theory and observations for normal-composition white dwarfs. For instance, the luminosity function of massive white dwarfs in the range of 0.9 < M_⋆/M_⊙ < 1.1 is relatively well fitted using current evolution models (Tremblay et al. 2019; Blouin et al. 2020). An additional multi-gigayear cooling delay for those objects can be ruled out. If a standard X(²²Ne) = 0.014 abundance is assumed instead of X(²²Ne) = 0.035, a ²²Ne-rich shell is formed, as in scenario (c) of Figure 3. Much less gravitational energy is then released, which, as needed, leads to shorter cooling delays (Table 1).

Instead of occurring at the onset of crystallization, those delays now occur once ∼60% of the core is already crystallized (assuming an homogeneous X(O) = 0.60 C/O initial profile). Very promisingly, a cooling delay of the order of ∼1 Gyr that manifests itself once ∼60% of the core is crystallized is missing from current C/O cooling tracks (Figure 2 of Blouin et al. 2020 and Figure 19 of Kilic et al. 2020). ²²Ne phase separation is very likely the solution to this second problem, although detailed evolutionary calculations that include ²²Ne gravitational settling in the liquid phase will be needed to confirm this hypothesis.

Finally, what about ²²Ne-rich normal-mass white dwarfs? A 0.6 M_⊙ white dwarf with X(²²Ne) = 0.035 and X(O) = 0.60 would experience a very long cooling delay of ≈12 Gyr (ΔB = 2.1 × 10⁴⁷ erg) during the formation of its ²²Ne-rich central core at $\mathrm{log}L/{L}_{\odot }\approx -3.9$. However, no discrepancy between the dynamical and photometric ages of the sort discovered for ultramassive white dwarfs (Cheng et al. 2019) has so far been identified for normal-mass objects. If confirmed, this nondetection can be explained in two ways:

• 1.  
  Very few normal-mass white dwarfs are enriched in ²²Ne to the level required to form a ²²Ne-rich central core. This would suggest that a negligible number of white dwarfs in the solar neighborhood were formed in α-rich environments, challenging Bauer et al.'s hypothesis that the delayed ultramassive population originates from such environments. In this context, the idea that the delayed ultramassive white dwarfs acquired their high ²²Ne content during a merger event (Blouin et al. 2020; Camisassa et al. 2020) may be more promising.
• 2.  
  Alternatively, the central layers of normal-mass white dwarfs may be much more rich in O than currently assumed (e.g., X(O) ≃ 0.85; Giammichele et al. 2018), which would lead to the formation of a ²²Ne-rich shell (Figure 3), with a much smaller cooling delay.

We have shown that ultramassive white dwarfs experiencing a long delayed cooling on the Q branch likely correspond to a population of ²²Ne-rich C/O white dwarfs in which ²²Ne phase separation induces a multi-gigayear cooling delay through the formation of a ²²Ne-rich central core. Such multi-gigayear delays are not observed for the vast majority of white dwarfs since, for canonical element abundances, ²²Ne phase separation leads instead to the formation of a ²²Ne-rich shell closer to the surface. In this case, a smaller cooling delay, which is actually required by current observational data, is generated. The fact that the qualitative outcome of ²²Ne phase separation depends on the ²²Ne mass fraction is the key attribute that allows this process to cause both small cooling delays at standard ²²Ne abundances and very long cooling delays at moderately above-average ²²Ne abundances. Detailed cooling sequences and population synthesis simulations will be needed to confirm our conclusions.

Future work should focus on building detailed evolutionary models to study the interplay of C/O phase separation and ²²Ne settling in the liquid phase with ²²Ne phase separation. The phase separation of ⁵⁶Fe, another minor species with a potentially large impact on white dwarf cooling (Xu & van Horn 1992), should also be investigated using modern simulation methods. Finally, this work motivates further efforts to better constrain the still poorly known C/O abundance profiles of white dwarfs (Giammichele et al. 2018), whose shape strongly affects the amount of energy released by ²²Ne phase separation.

We are grateful to the anonymous referee for useful comments that have improved the manuscript. S.B. thanks A. Bédard and P. Brassard for useful discussions on the STELUM code. Research presented in this article was supported by the Laboratory Directed Research and Development program of Los Alamos National Laboratory under project number 20190624PRD2. This work was performed under the auspices of the U.S. Department of Energy under Contract No. 89233218CNA000001.