Buoyant crystals halt the cooling of white dwarf stars – Nature

Cooling models: basic assumptions

Our C–O white dwarf cooling models were computed using the STELUM code²³. We generated and evolved 1.00–1.25 M_☉ model structures consisting of a ¹²C–¹⁶O–²²Ne core surrounded by a He mantle and an outermost H layer. We adopted envelope layer masses M(H) = 10⁻¹⁰ M_★ and M(He) = 10⁻⁶ M_★, along with an initially uniform core with mass-fraction abundances X(¹⁶O) = 0.50 and X(²²Ne) = 0.03. From now on, we refer to this as our post-merger chemical structure. The envelope composition is based on empirical characterization of Q-branch white dwarfs^9,29 and is also consistent with the expected outcome of nuclear burning following a stellar merger³⁰. The core composition is taken from recent theoretical work on white dwarf–subgiant mergers¹⁰, with the exception that ²²Ne is used as a proxy for all neutron-rich species. This choice is motivated by the fact that white dwarf–subgiant mergers represent the most promising scenario for the formation of ultramassive white dwarfs with both a C–O core and a high content of neutron-rich impurities. The uniform composition is a direct consequence of the merger event, which generates thermohaline mixing that homogenizes the entire core¹⁰. This also explains the slightly lower central ¹⁶O abundance compared with standard values expected from single-star evolution^16,37, X(¹⁶O) ≃ 0.55–0.60, as the ¹⁶O normally concentrated towards the centre is homogeneously redistributed over the core.

We ignored general relativistic effects, which is an excellent approximation in the stellar mass range considered^38,39. We also omitted rotation and magnetic fields; although some merger remnants are observed to be rapidly rotating and/or strongly magnetic^40,41, the effects on the stellar structure and evolution are expected to be negligible^2,42. We used the radiative opacities of the OPAL project⁴³ and the latest conductive opacities^44,45. Convective energy transport in the envelope was handled following the ML2 version of the mixing-length theory^46,47. We used a simple grey atmosphere as the outer boundary condition, a valid approximation for Q-branch white dwarfs as they have not yet undergone surface–core convective coupling^12,48.

Element transport was modelled through the usual scheme in which the time-dependent transport equations are fully coupled to the stellar-structure equations²³. Standard microscopic diffusion was included in the gas and liquid phases using newly implemented diffusion coefficients^49,50. This allows for a state-of-the-art description of the heating effect owing to ²²Ne diffusion in the core^18,26,51. We considered gravitational settling and chemical diffusion but ignored thermal diffusion, as it is negligible in the nearly isothermal core of white dwarfs. We omitted chemical mixing resulting from convection in the envelope because this process can give rise to severe numerical difficulties⁵². This assumption is certainly valid for the DA-type members of the Q-branch population, whose H-dominated atmosphere indicates that the superficial H layer has not been altered by convection. This is, however, not true for their DQ-type counterparts, which probably owe their H-deficient surface composition to the convective mixing of the outer H, He and C layers^29,52. Nevertheless, we verified that artificially imposing a uniform He–C envelope with typical abundances^29,53 barely affects the evolution, whether before, during or after the distillation phase. The overall opacity turns out to be similar to that provided by a thin pure-H layer, which thus represents an acceptable approximation even for the DQ stars.

Cooling models: crystallization and distillation

As a prerequisite to the inclusion of the distillation mechanism in STELUM, we first upgraded the treatment of standard C–O fractionation during crystallization in the absence of neutron-rich impurities. We replaced the approximate heating term used previously^23,54 with a detailed description of the composition changes in the liquid and solid regions, similar to that adopted in the MESA code¹⁶. This procedure involves the following steps, which are performed at the end of each evolutionary time step:

1. 1.

   Identifying the layers that have solidified during the last time step. The phase transition is assumed to occur when the plasma coupling parameter Γ becomes larger than the critical value Γ_cr predicted by the C–O phase diagram of ref. ⁵⁵. Note that the coupling parameter of the plasma mixture is defined as Γ = ⟨Z^5/3⟩e²/a_ek_BT, in which ⟨Z^5/3⟩ = ∑_i Z_i^5/3x_i (x_i being the number fraction of the ionic species of charge Z_i), e is the elementary charge, a_e is the radius of a sphere whose volume corresponds to the mean volume per electron, k_B is the Boltzmann constant and T is the temperature.

2. 2.

   Adjusting the composition of the newly crystallized layers. The ¹⁶O abundance is increased following the same phase diagram as above, whereas the ¹²C abundance is decreased by an equal amount. The composition of the solid region is then held fixed for the rest of the evolutionary calculation.

3. 3.

   Adjusting the composition of the innermost liquid layer to enforce total elemental mass conservation. Given the abundance changes in the solid region, the first liquid layer above the crystallization front is enriched in ¹²C and depleted in ¹⁶O, which produces an inverted molecular-weight gradient.

4. 4.

   Imposing a uniform composition in the dynamically unstable liquid layers. The instability arising from the molecular-weight gradient is assumed to mix material outward in a homogeneous and instantaneous way^16,56,57,58. The extent of the mixed region is found by redistributing the elements over an increasing number of layers until the Ledoux stability criterion is satisfied^16,58,59.

In the case of ²²Ne distillation, fractionation produces ²²Ne-poor crystals that float up and eventually melt, thereby displacing ²²Ne-rich liquid downward. The net result is that the central crystal-forming region remains globally liquid (hence the term crystal-forming rather than crystallized) and gradually becomes enriched in ¹²C and ²²Ne and depleted in ¹⁶O at the expense of the outer layers⁸. We implemented this process in STELUM using the same algorithm as above but with a few modifications to steps 1 and 2. In step 1, the impact of ²²Ne on the onset of crystallization is taken into account through a new analytical fit to the high-resolution C–O–Ne phase diagram of ref. ⁸. The fitting formula specifies the critical coupling parameter of C, Γ_cr,C, and is expressed as a correction to the C–O phase diagram of ref. ⁵⁵,

in which Γ_cr,C⁰ is the uncorrected value, x_O and x_Ne are the number fractions of ¹⁶O and ²²Ne, respectively, and the numerical coefficients are c₁ = 1,096.69, c₂ = −3,410.33 and c₃ = 2,408.44. Because the C–O–Ne phase diagram of ref. ⁸ is limited to relatively small ²²Ne abundances, the above expression is only valid for x_Ne < 0.04. Furthermore, the crystal-forming layers so identified are assumed to remain liquid, meaning that their composition is allowed to change at each subsequent time step. In step 2, the chemical evolution of the crystal-forming region is handled using a prescription for the elemental abundances as a function of the coupling parameter of C, Γ_C. Current phase diagrams specify the composition at the beginning and end of distillation^8,24, but the exact trajectory followed between these two points remains unknown. We opted for the simplest possible prescription, in which the abundances change linearly with Γ_C between the initial and final states. The latter is given by (Γ_C, x_C, x_O, x_Ne) = (208, 0.8, 0.0, 0.2) in terms of number fractions, corresponding to a ²²Ne mass fraction X(²²Ne) ≃ 0.31 (ref. ⁸). Once the composition of the crystal-forming layers has been updated, steps 3 and 4 are applied as described above, except that ²²Ne mass conservation and redistribution are also considered. The energy released by ²²Ne distillation is naturally taken into account through the coupling of the chemical-transport and stellar-structure equations in the STELUM code²³.

The distillation process is assumed to stop once the whole ²²Ne content of the star has been transported to the central crystal-forming region. This is admittedly an uncertain aspect of our modelling, as it is entirely possible that distillation terminates earlier; this point is addressed in the ‘Sensitivity to the distillation implementation’ section. At this stage, the outer ²²Ne-free layers of the core are expected to start crystallizing normally, as there is no more ²²Ne to produce a density reversal with respect to the overlying liquid. These stable solid layers presumably block the floating motion of the buoyant crystals immediately underneath, which thus remain in place and in turn obstruct crystals originating from deeper layers, and so on. The probable outcome is that the whole ²²Ne-rich region solidifies at once, thereby completely halting the distillation process. For this reason, in our models, the inner ²²Ne-rich core is assumed to freeze in place at the end of distillation, and its composition is held fixed afterwards. This region typically encompasses ≃25% of the total mass of the star; the sudden solidification of that much mass releases a large amount of latent heat, which can cause notable numerical instabilities. To avoid such instabilities, the progression of the crystallization front is artificially slowed down by setting the critical coupling parameter of the entire core to its value at the centre (where it is highest, as the ²²Ne abundance is highest^24,60). This ensures that the energy release is more gradual in both space and time, thereby greatly improving numerical convergence. This procedure slightly accelerates the post-distillation cooling of our models, but this is of no consequence to our analysis given that most delayed white dwarfs are not old enough to have reached this phase (Fig. 4a). Once the crystallization front reaches the outer ²²Ne-free layers, standard C–O fractionation occurs as described above.

To better illustrate our modelling assumptions, we show in Extended Data Fig. 1 the composition profile of our 1.15 M_☉ model just before, during and just after distillation, as well as at the very end of the evolutionary sequence (at which point nearly the entire core has solidified). The distillation mechanism produces an increasingly ²²Ne-rich, ¹⁶O-poor inner core, which also expands with time as more layers start forming buoyant crystals. In this region, the abundance profiles simply reflect the Γ_C profile of the star² as a result of our chemical-evolution prescription. At the end of distillation, there is no more ²²Ne in the outer portion of the core and the composition profile of the inner ²²Ne-rich layers (m/M_★ ≲ 0.25) is subsequently held fixed. In the ²²Ne-free region, the final ¹⁶O abundance profile is the well-known result of standard fractionation during crystallization^14,16.

During the distillation process, the gravitational energy released almost exactly balances the energy radiated away, such that the surface luminosity and temperature remain essentially constant. In practice, the use of a finite time step in our numerical calculations gives rise to small oscillations around the mean values (typically about 1% in T_eff). We mitigated these numerical artefacts by applying a simple moving-average smoothing procedure (with an averaging window of around 1 Gyr) along the luminosity–temperature plateau of our evolutionary sequences. This is important to ensure an accurate and well-behaved interpolation of the cooling models in the population-synthesis simulations. All cooling sequences used and shown in this paper are the final smoothed versions.

A key characteristic of the observed Q branch reproduced by our calculations is the fact that it is much narrower than the luminosity range of crystallization itself⁵. This is explained by two properties of the distillation mechanism. First, distillation operates only in the initial phase of the crystallization process, as it must stop once all neutron-rich impurities have been displaced to the central layers. Second, the rate of energy release during the distillation process itself is not uniform. It is roughly proportional to the total mass of the coexistence region and therefore peaks only after distillation has operated for some time. This explains the offset between the start of distillation and the cooling pause in Fig. 4a.

Population synthesis

We considered three distinct ultramassive white dwarf populations in our population synthesis: 5–9% of C–O white dwarfs that undergo distillation, ≃30% of O–Ne white dwarfs with merger delays and the remaining standard O–Ne white dwarfs. The first group is characterized by a composition similar to that of white dwarf–subgiant merger remnants¹⁰, although we do not exclude other evolutionary channels. We assigned a proportion of 5–9% to this population as constrained by a previous kinematic analysis⁵. These objects were modelled using the cooling calculations detailed above. The second group consists of double white dwarf merger products. We supposed that it accounts for ≃30% of the ultramassive population (with a small dependence on the mass)⁶¹ and that they harbour O–Ne cores⁶². Finally, the remainder are standard O–Ne white dwarfs resulting from single-star evolution and were modelled using existing cooling tracks²¹.

We assumed that the age distribution of stars in the 150-pc solar neighbourhood is constant between 0 and 10.5 Gyr (ref. ³⁵). We did not explicitly model the initial mass function, because a sizable fraction of the ultramassive population emerged from stellar mergers. Instead, we used a large sample of white dwarfs¹⁷ to infer a ∝M_★⁻² empirical probability distribution function for white dwarfs more massive than 1 M_☉. For a given white dwarf mass taken from this distribution, we then added the pre-white dwarf lifetime^63,64 and a merger delay when applicable. We modelled the delay-time distribution using a log-normal distribution. For the natural logarithm of ages expressed in Gyr, μ = 0.5 and σ = 0.7, which corresponds to an average delay of 2.1 Gyr (ref. ⁶¹). The Gaia magnitudes were then calculated using pure-H-atmosphere models^65,66, except for 50% of the extra-delayed C–O core population, for which we instead used He-dominated, C-polluted atmosphere models with N(C)/N(He) = 0.1 (ref. ⁵³). This assumption is motivated by a kinematic analysis of spectroscopically observed Q-branch white dwarfs, which revealed that the extra-delayed population is composed of DA and DQ stars in similar proportions, whereas the remainder comprises very few (if any) DQ stars⁵. For consistency with previous work⁵, we only considered ultramassive white dwarfs whose Gaia photometry indicates a mass comprised in the 1.08–1.23 M_☉ interval when modelled assuming H atmospheres and O–Ne cores.

Model atmospheres

An existing DQ model atmosphere grid⁵³ was extended to calculate the Gaia magnitudes of Q-branch objects with C-polluted photospheres. The original model grid stopped at T_eff = 16,000 K and log g = 9.0; we expanded it to T_eff = 20,000 K and log g = 9.5 using the same model atmosphere code^67,68.

Sample selection

The white dwarf sample shown in Figs. 1 and 3 was taken from the Gaia EDR3 white dwarf catalogue⁶⁹. We selected all high-confidence white dwarf candidates (that is, those with a probability of being white dwarfs of at least 90%, P_WD ≥ 0.9) that lie within D = 150 pc of the Sun. This distance cut-off is justified in Extended Data Fig. 2, in which we show the cumulative distribution functions (CDFs) of the distances of high-confidence white dwarf candidates for three different M_G bins. For the faintest bin relevant to our analysis (13.5 < M_G ≤ 14.5), the CDF departs from the expected ∝D³ relation at ≃150 pc. This indicates that the sample starts to be noticeably incomplete at faint magnitudes past 150 pc. Accordingly, we only selected stars within 150 pc to avoid biasing the sample in favour of brighter stars.

We cross-matched all 458 white dwarfs from the 150-pc sample that are used in our analysis (that is, those falling between the dashed grey lines in Fig. 1) with the Montreal White Dwarf Database⁷⁰. We found that only three out of the 79 objects with a known spectral type are of the DB class, with DAs and DQs making up most of the remainder. This justifies our use of pure-H-atmosphere models for all objects in our simulations, except for half of the extra-delayed C–O core population⁵, which we modelled using DQ model atmospheres.

Sensitivity to the assumed composition

As mentioned above, our nominal cooling sequences assume a post-merger chemical structure, with an extremely thin envelope and an initially uniform core. To investigate the impact of these assumptions on our results, we also computed a set of white dwarf models with more standard envelope and core stratifications (that is, emulating the predictions of single-star-evolution calculations). For this test, we used M(H) = 10⁻⁶ M_★, M(He) = 10⁻³ M_★ (refs. ^21,28) and the non-uniform C–O profile of the 1.0 M_☉ model of ref. ¹⁶ (but still assuming X(²²Ne) = 0.03). The resulting cooling sequences (Extended Data Fig. 3a) are substantially different from our nominal case: for a given mass, the constant-luminosity phase owing to ²²Ne distillation occurs at a higher luminosity and is thus shorter. The reason is that the models start crystallizing earlier, a behaviour owing mainly to the larger He layer mass, but also to the higher central O abundance. In particular, the larger amount of He increases the overall transparency of the envelope, resulting in a higher surface luminosity for a given core temperature. The earlier and shorter cooling delays lead to a much poorer agreement between the simulated and observed Q-branch overdensities, in terms of both location and amplitude (Extended Data Fig. 4a). The thin He layer of our post-merger-type models also explains the longer distillation cooling delays compared with previous estimates⁸.

Another potentially important parameter for our study is the initial ²²Ne abundance, which dictates the amount of gravitational energy available for release through distillation. Given current uncertainties on this quantity¹⁰, we computed further cooling sequences assuming a lower ²²Ne content of X(²²Ne) = 0.025 along with our default post-merger chemical structure. For a given mass, the time evolution of the surface luminosity is qualitatively similar to the nominal case, but the cooling delay owing to distillation is, unsurprisingly, smaller (≃6–11 Gyr for 1.00–1.25 M_☉). Repeating our population-synthesis simulation with these models, we find that the results are only weakly affected (Extended Data Fig. 4b). Similarly, we found that increasing the ²²Ne content to X(²²Ne) = 0.04 has a negligible impact on the predicted Q branch. The cooling delays for X(²²Ne) = 0.03 are already long enough that most delayed white dwarfs are still stuck on the Q branch at the final 10.5-Gyr age.

Sensitivity to the distillation implementation

The C–O–Ne phase diagram predicts that the solid phase is buoyant as long as the abundance of ²²Ne in the liquid remains above a well-defined threshold⁸. As the distillation process unfolds, there is inevitably a point at which the ²²Ne abundance in the outermost crystal-forming layer reaches this critical value. What happens past that juncture remains unclear. A priori, distillation could stop because the solid phase is no longer buoyant. However, despite not being buoyant, the solid phase is still depleted in ²²Ne. This means that crystallization should increase the ²²Ne abundance in the liquid, thereby reinstating the conditions required for distillation. Stopping distillation would then be an unstable state. To address this uncertainty, we calculated further cooling sequences corresponding to an extreme scenario in which distillation is terminated when the threshold abundance is first reached (Extended Data Fig. 3b). The distillation-induced cooling delay is then shorter, but it remains sufficient to explain the Q-branch overdensity (Extended Data Fig. 5). Our conclusions are therefore insensitive to the treatment of distillation past this critical juncture.

Microscopic diffusion of ²²Ne

It has been claimed that simple microscopic diffusion (often referred to as gravitational settling) of ²²Ne in C–O white dwarfs can possibly explain the Q-branch cooling anomaly¹⁹, an idea similar to the ‘sedimentars’ proposed two decades earlier⁷¹. However, this mechanism cannot reproduce the narrowness of the Q branch (see Fig. 4 of ref. ¹⁹). Moreover, to obtain a long enough delay time, a very large ²²Ne mass fraction of 0.06 is required, a value that far surpasses what can be explained by any known evolutionary pathway. This scenario also requires that as much as about 50% of all ultramassive white dwarfs experience a delayed cooling, in sharp disagreement with kinematic constraints⁵. In short, this is not a viable scenario.

Our cooling models are consistent with these conclusions. Ignoring distillation and assuming X(²²Ne) = 0.06 (along with our usual post-merger composition), we find that the slowdown of the cooling process is much less pronounced than in our default setup (Extended Data Fig. 3c). Gravitational settling is simply far less efficient than distillation at transporting neutron-rich species to the centre of the star (Extended Data Fig. 6). We note that the settling-induced cooling delays of our X(²²Ne) = 0.06 models are slightly shorter than those reported by ref. ¹⁹, in line with independent calculations¹⁶. A population-synthesis simulation in which 5–9% of ultramassive white dwarfs follow this evolutionary scenario yields a small bump centred on the Q branch, but this overdensity is too shallow and wide to reproduce the observations (Extended Data Fig. 7). This is a direct consequence of the fact that, unlike distillation, gravitational settling operates over a wide range of luminosities (Extended Data Fig. 3c).

Sensitivity to star-formation history

We assumed a uniform stellar-age distribution for our 150-pc solar neighbourhood sample for up to 10.5 Gyr. This distribution is both consistent with the recent analysis of a nearby white dwarf sample³⁵ and expected because of the counterbalancing of the declining star-formation rate⁷² by the thinner disk scale height for younger stars⁷³. However, independent studies instead find strong variations in the effective star-formation history of the Galactic disk, with a stellar formation burst 2–3 Gyr ago⁷⁴. To verify whether this could have an effect on our conclusions²⁰, we show in Extended Data Fig. 8 the results of population-synthesis simulations identical to those presented in Fig. 3 but using a non-uniform age distribution inspired by the results of ref. ⁷⁴. The age distribution is assumed to correspond to the sum of a uniform distribution and a Gaussian function centred at 2.5 Gyr ago with σ = 1 Gyr and an amplitude that increases the number of 2.5-Gyr-old stars by 60% compared with the uniform baseline. A small bump associated with the stellar-formation burst appears just after the Q-branch peak, but its amplitude and location do not match the Q-branch signal. The agreement between our nominal simulation and the Gaia data has worsened compared with Fig. 3, but it is clear that distillation is still needed to explain the sharp peak between ζ = 13.0 and 13.2. Our conclusions are therefore insensitive to the assumed stellar formation history.

Sensitivity to merger delay

The log-normal distribution for merger time delays used in our population synthesis mimics the results of ref. ⁶¹ for double white dwarf mergers. However, this distribution may not be appropriate for all types of merger that can contribute to the ultramassive population. We repeated the population-synthesis simulations of Fig. 3 assuming no merger delay for all stars and found that our results are barely affected (Extended Data Fig. 9a). Another test was performed in which the delay times were doubled for all merger remnants, and the same conclusion was reached (Extended Data Fig. 9b). Our results are therefore not sensitive to the assumed merger delay time distribution. This is expected because a merger creates only an overall shift of cooling times but no pile-up along the cooling track.

Population density at ζ > 13.4

A tension between the observed and simulated populations appears at ζ > 13.4 in Fig. 3. In the 13.4 < ζ < 13.8 range, this may be because of the high conductive opacities assumed in the O–Ne white dwarf models²¹ used for the bulk of the synthetic population^45,75. At ζ > 13.8, the observed population is probably contaminated by faint, lower-mass white dwarfs with uncertain Gaia data (see the bottom-right corner of Fig. 1). In any case, this tension does not affect our conclusions about the Q branch and distillation, which are based on the analysis of white dwarfs with ζ < 13.4.

Kinematics of Q-branch white dwarfs

It has been claimed that the velocity dispersion of white dwarfs on the Q branch in the direction of Galactic rotation is higher than even the oldest population of the Galactic disk⁷⁶. However, we believe that is only an apparent mismatch owing to the invalidity of fitting a Gaussian tail to a non-Gaussian velocity distribution when its dispersion is large. Below we show that the Q-branch kinematics matches that of stars in the solar vicinity with ages older than 1–2 Gyr, which is a prediction from the cooling-delay scenario.

To avoid inconsistencies caused by selection effects in different datasets used in the literature, we directly compared the velocity distribution of Q-branch white dwarfs to a sample of M dwarfs in Gaia, with the same distance cut of 150 pc from the Sun and a colour cut of 1.49 ≤ BP − RP ≤ 1.50. The absolute magnitude of this sample is higher than the Q branch (M_G ≃ 7.5) and thus it is also complete at the same distance. In a cooling-delay scenario, the Q branch is expected to have the same age distribution as M dwarfs older than 1–2 Gyr, because both samples have no age bias or selection in that range. We computed the velocity dispersion of both samples projected in the Galactic radial, rotation and vertical directions, after correcting for the contribution from young objects. We first removed low-velocity (<20 km s⁻¹) objects, that is, young stars from both samples. The M-dwarf sample then has dispersions of 34, 31 and 25 km s⁻¹, respectively, whereas the Q branch (ζ within 13.0–13.2, mass within 1.08–1.23 M_☉) has dispersions of 31, 30 and 23 km s⁻¹. The Q-branch dispersion is slightly smaller than the M-dwarf sample, but because the Q-branch sample contains a high fraction of normal white dwarfs (about half, as estimated in ref. ⁵ or read from Fig. 3) that are 1–2 Gyr old, the tail of their velocity distribution may still contaminate the range above 20 km s⁻¹. We estimated their contribution using the velocities of white dwarfs with 12.4 < ζ < 12.8. After such a correction, the velocity dispersion of the Q branch becomes 32, 32 and 25 km s⁻¹, accurately matching the M-dwarf sample and supporting the idea that the Q branch is an accumulation of stars with a wide range of ages caused by a long cooling delay. It may be worth further investigating the slightly higher ratio between dispersions along the Galactic rotation and radial directions, but we find no evidence that the Q-branch velocity dispersion exceeds the expected disk values.