Strange IndiaStrange India

Cooling models: basic assumptions

Our C–O white dwarf cooling models were computed using the STELUM code23. We generated and evolved 1.00–1.25 M model structures consisting of a 12C–16O–22Ne core surrounded by a He mantle and an outermost H layer. We adopted envelope layer masses M(H) = 10−10M and M(He) = 10−6M, along with an initially uniform core with mass-fraction abundances X(16O) = 0.50 and X(22Ne) = 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 dwarfs9,29 and is also consistent with the expected outcome of nuclear burning following a stellar merger30. The core composition is taken from recent theoretical work on white dwarf–subgiant mergers10, with the exception that 22Ne 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 core10. This also explains the slightly lower central 16O abundance compared with standard values expected from single-star evolution16,37, X(16O)  0.55–0.60, as the 16O 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 considered38,39. We also omitted rotation and magnetic fields; although some merger remnants are observed to be rapidly rotating and/or strongly magnetic40,41, the effects on the stellar structure and evolution are expected to be negligible2,42. We used the radiative opacities of the OPAL project43 and the latest conductive opacities44,45. Convective energy transport in the envelope was handled following the ML2 version of the mixing-length theory46,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 coupling12,48.

Element transport was modelled through the usual scheme in which the time-dependent transport equations are fully coupled to the stellar-structure equations23. Standard microscopic diffusion was included in the gas and liquid phases using newly implemented diffusion coefficients49,50. This allows for a state-of-the-art description of the heating effect owing to 22Ne diffusion in the core18,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 difficulties52. 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 layers29,52. Nevertheless, we verified that artificially imposing a uniform He–C envelope with typical abundances29,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 previously23,54 with a detailed description of the composition changes in the liquid and solid regions, similar to that adopted in the MESA code16. 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. 55. Note that the coupling parameter of the plasma mixture is defined as Γ = Z5/3e2/aekBT, in which Z5/3 = ∑i Zi5/3xi (xi being the number fraction of the ionic species of charge Zi), e is the elementary charge, ae is the radius of a sphere whose volume corresponds to the mean volume per electron, kB is the Boltzmann constant and T is the temperature.

  2. 2.

    Adjusting the composition of the newly crystallized layers. The 16O abundance is increased following the same phase diagram as above, whereas the 12C 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 12C and depleted in 16O, 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 way16,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 satisfied16,58,59.

In the case of 22Ne distillation, fractionation produces 22Ne-poor crystals that float up and eventually melt, thereby displacing 22Ne-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 12C and 22Ne and depleted in 16O at the expense of the outer layers8. 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 22Ne 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. 8. 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. 55,

$${\Gamma }_{{\rm{cr}},{\rm{C}}}={{\Gamma }_{{\rm{cr}},{\rm{C}}}}^{0}+\left({c}_{1}{x}_{{\rm{O}}}+{c}_{2}{{x}_{{\rm{O}}}}^{2}+{c}_{3}{{x}_{{\rm{O}}}}^{3}\right){x}_{{\rm{Ne}}},$$

in which Γcr,C0 is the uncorrected value, xO and xNe are the number fractions of 16O and 22Ne, respectively, and the numerical coefficients are c1 = 1,096.69, c2 = −3,410.33 and c3 = 2,408.44. Because the C–O–Ne phase diagram of ref. 8 is limited to relatively small 22Ne abundances, the above expression is only valid for xNe < 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 distillation8,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, xC, xO, xNe) = (208, 0.8, 0.0, 0.2) in terms of number fractions, corresponding to a 22Ne mass fraction X(22Ne)  0.31 (ref. 8). Once the composition of the crystal-forming layers has been updated, steps 3 and 4 are applied as described above, except that 22Ne mass conservation and redistribution are also considered. The energy released by 22Ne distillation is naturally taken into account through the coupling of the chemical-transport and stellar-structure equations in the STELUM code23.

The distillation process is assumed to stop once the whole 22Ne 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 22Ne-free layers of the core are expected to start crystallizing normally, as there is no more 22Ne 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 22Ne-rich region solidifies at once, thereby completely halting the distillation process. For this reason, in our models, the inner 22Ne-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 22Ne abundance is highest24,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 22Ne-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 22Ne-rich, 16O-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 star2 as a result of our chemical-evolution prescription. At the end of distillation, there is no more 22Ne in the outer portion of the core and the composition profile of the inner 22Ne-rich layers (m/M 0.25) is subsequently held fixed. In the 22Ne-free region, the final 16O abundance profile is the well-known result of standard fractionation during crystallization14,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 Teff). 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 itself5. 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 remnants10, although we do not exclude other evolutionary channels. We assigned a proportion of 5–9% to this population as constrained by a previous kinematic analysis5. 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)61 and that they harbour O–Ne cores62. Finally, the remainder are standard O–Ne white dwarfs resulting from single-star evolution and were modelled using existing cooling tracks21.

We assumed that the age distribution of stars in the 150-pc solar neighbourhood is constant between 0 and 10.5 Gyr (ref. 35). 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 dwarfs17 to infer a M−2 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 lifetime63,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. 61). The Gaia magnitudes were then calculated using pure-H-atmosphere models65,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. 53). 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 stars5. For consistency with previous work5, 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 grid53 was extended to calculate the Gaia magnitudes of Q-branch objects with C-polluted photospheres. The original model grid stopped at Teff = 16,000 K and log g = 9.0; we expanded it to Teff = 20,000 K and log g = 9.5 using the same model atmosphere code67,68.

Sample selection

The white dwarf sample shown in Figs. 1 and 3 was taken from the Gaia EDR3 white dwarf catalogue69. We selected all high-confidence white dwarf candidates (that is, those with a probability of being white dwarfs of at least 90%, PWD ≥ 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 MG bins. For the faintest bin relevant to our analysis (13.5 < MG ≤ 14.5), the CDF departs from the expected D3 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 Database70. 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 population5, 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−6M, M(He) = 10−3M (refs. 21,28) and the non-uniform C–O profile of the 1.0 M model of ref. 16 (but still assuming X(22Ne) = 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 22Ne 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 estimates8.

Another potentially important parameter for our study is the initial 22Ne abundance, which dictates the amount of gravitational energy available for release through distillation. Given current uncertainties on this quantity10, we computed further cooling sequences assuming a lower 22Ne content of X(22Ne) = 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 22Ne content to X(22Ne) = 0.04 has a negligible impact on the predicted Q branch. The cooling delays for X(22Ne) = 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 22Ne in the liquid remains above a well-defined threshold8. As the distillation process unfolds, there is inevitably a point at which the 22Ne 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 22Ne. This means that crystallization should increase the 22Ne 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 22Ne

It has been claimed that simple microscopic diffusion (often referred to as gravitational settling) of 22Ne in C–O white dwarfs can possibly explain the Q-branch cooling anomaly19, an idea similar to the ‘sedimentars’ proposed two decades earlier71. However, this mechanism cannot reproduce the narrowness of the Q branch (see Fig. 4 of ref. 19). Moreover, to obtain a long enough delay time, a very large 22Ne 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 constraints5. In short, this is not a viable scenario.

Our cooling models are consistent with these conclusions. Ignoring distillation and assuming X(22Ne) = 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(22Ne) = 0.06 models are slightly shorter than those reported by ref. 19, in line with independent calculations16. 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 sample35 and expected because of the counterbalancing of the declining star-formation rate72 by the thinner disk scale height for younger stars73. 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 ago74. To verify whether this could have an effect on our conclusions20, 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. 74. 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. 61 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 models21 used for the bulk of the synthetic population45,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 disk76. 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 (MG 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−1) objects, that is, young stars from both samples. The M-dwarf sample then has dispersions of 34, 31 and 25 km s−1, 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−1. 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. 5 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−1. 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−1, 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.

Source link


Leave a Reply

Your email address will not be published. Required fields are marked *