© The Authors 2024

1. Introduction

White dwarfs (WDs) with either carbon-oxygen (CO) or oxygen-neon (ONe) cores are the most common end-stage of the evolution of single stars with initial masses of up to ∼8–9 M_⊙. Given the current age of the universe and the shape of the stellar initial mass function, the overwhelming majority of existing WDs produced by single stars (and non-interacting binaries) are the progeny of objects with initial masses typically between ∼0.8–1.0 and ∼6–7 M_⊙ that have a CO core.

The evolution of WDs is a cooling process (see e.g. the review by Saumon et al. 2022) that produces a well-defined relationship between a WD’s luminosity and its cooling age, which can be employed as a useful cosmic clock. In the last 2–3 decades, theoretical models of CO-core WDs, together with photometric, spectroscopic, and asteroseismic data, have been employed to constrain the star formation history of the local disk from its WD population (e.g. Winget et al. 1987; Oswalt et al. 1996; Torres & García-Berro 2016; Kilic et al. 2017; Tononi et al. 2019; Cukanovaite et al. 2023), the ages of several open clusters (e.g. Richer et al. 1998; von Hippel 2005; Bedin et al. 2008a, 2010, 2015; Bellini et al. 2010; García-Berro et al. 2010), and globular clusters (e.g. Hansen et al. 2004, 2007; Winget et al. 2009; Bedin et al. 2009, 2019; Goldsbury et al. 2012). They have also been used as proxies for laboratories to investigate open questions in theoretical physics (see e.g. the review by Isern et al. 2022, and references therein).

The accuracy of the ages derived from WD cosmo-chronology rests on the accuracy of WD models, and hence on the correct description of all the processes that can contribute to a WD energy reservoir. There has been much interest in processes beyond the crystallization and phase separation of the CO chemical mixture (see e.g. Segretain et al. 1994; Blouin & Daligault 2021, and references therein) that can contribute substantially to the CO WD energy budget – and hence cooling times – and involve the presence of ²²Ne, the most abundant metal in the core of a CO WD after ¹²C and ¹⁶O (all other metals are at least about one order of magnitude less abundant than ²²Ne in the CO cores). This neon isotope is produced during the He-burning evolution of the WD progenitors, and its abundance across the CO core is to a very good approximation constant, roughly equal (in mass fraction) to the progenitors’ initial total metallicity (see e.g. Salaris et al. 2022, and references therein), typically about 1.5–2% for an initial solar chemical composition.

Given that ²²Ne nuclei have a larger mass-to-charge ratio than the dominant ¹²C and ¹⁶O components, this results in a downward gravitational force on ²²Ne and a slow diffusion towards the centre in the liquid layers, with the release of gravitational energy (see Bildsten & Hall 2001; Deloye & Bildsten 2002; García-Berro et al. 2008; Camisassa et al. 2016; Salaris et al. 2022; Bauer 2023). This extra energy contribution can appreciably impact the models’ cooling times and depends on the progenitor initial metallicity: higher initial metallicities imply larger ²²Ne mass fractions and longer cooling time delays due to diffusion.

Very recently Blouin et al. (2021) performed state-of-the-art Monte Carlo plasma simulations to study the phase separation of ²²Ne in a crystallizing CO WD. Based on their results, they revived the idea of ²²Ne distillation (Isern et al. 1991; Segretain 1996), where ²²Ne is efficiently transported towards the WD centre and a potentially substantial amount of gravitational energy is released. This transport mechanism is triggered by the partial exclusion of ²²Ne from the solid phase. Under the right conditions, this can render the solid buoyant by making it less dense than the coexisting liquid. The buoyant crystals then float up (and melt), displacing heavier liquid towards the centre of the star. This macroscopic transport of ²²Ne can be much more efficient than ²²Ne diffusion in the liquid phase, leading to larger effects on WD evolution and cooling delays of up to several Gyr.

While neon diffusion has been included in the computation of WD cooling models (see e.g. García-Berro et al. 2008; Althaus et al. 2010; Camisassa et al. 2016; Salaris et al. 2022; Bauer 2023), the recently studied distillation process has yet to be modelled within full evolutionary calculations of WDs. In this paper we present the first calculations that include CO crystallization and phase separation plus ²²Ne diffusion in the liquid phase and distillation during CO crystallization. We calculated WD models (and isochrones) for an initial progenitor metallicity (about twice solar) appropriate for the super-solar old open cluster NGC 6791, which has already been studied as a test bench for the effect of neon diffusion on its WD age determination (e.g. García-Berro et al. 2010). We compare our calculations with the cluster cooling sequence to assess whether models that include ²²Ne distillation are consistent with observations. This comparison constitutes a strong test for the models and the efficiency of neon distillation in crystallizing WDs, given the old age of the cluster (and hence the presence of WDs undergoing crystallization), its very high metallicity (and hence a high ²²Ne abundance, equal to about 3% by mass, which maximizes the effect of distillation), and the existing very tight constraints on the cluster age and distance.

In Sect. 2 we re-examine the cluster age determination from its WDs in light of the constraints set by the study of several of its eclipsing binaries and the more recent WD calculations for the effect of ²²Ne diffusion. Section 3 presents our new WD models, which include ²²Ne distillation, and highlights the impact of this process on the models’ chemical stratification and cooling times. A comparison of theoretical WD isochrones and luminosity functions (LFs) calculated from our new WD models with the observed NGC 6791 counterparts follows in Sect. 4. In Sect. 5 we propose and describe an independent observational test of the efficiency of neon distillation in this cluster’s WDs; this test, which would make use of the James Webb Space Telescope (JWST), could potentially help resolve the current uncertainties on the electron conduction opacities in the regime of WD envelopes. A brief summary and conclusions are provided in Sect. 6.

2. The cooling sequence of NGC 6791

NGC 6791 is one of the richest open clusters, is unusually old (age around 8 Gyr), and is very metal-rich ([Fe/H] ∼ +0.3). Its colour-magnitude diagram (CMD) strongly resembles that of a Galactic globular cluster.

By means of Hubble Space Telescope (HST) ACS imaging, Bedin et al. (2008a) reached the end of the cluster’s well-populated WD cooling sequence, whose differential LF (see Fig. 1) displays a peak at magnitude m_F606W = 27.45 ± 0.05 and a fainter peak at m_F606W = 28.15 ± 0.05. While the fainter peak is associated with the termination of the cooling sequence, and hence is the cluster age indicator, there is still no agreement on the origin of the brighter peak (see e.g. Hansen 2005; Bedin et al. 2008b, for two very different explanations).

Fig. 1. The cooling sequence of NGC 6791. Upper panel: observed WD LF of NGC 6791 (corrected for completeness) compared with a theoretical LF calculated for an age of 8.5 Gyr using a [Fe/H] = 0.3 WD isochrone (with c07 opacities) from models calculated without including the effect of 22Ne distillation (see the main text for details). Lower panel: CMD of the observed WD cooling sequence of NGC 6791 together with [Fe/H] = 0.3, 8.5 Gyr isochrones calculated without 22Ne distillation and with either the c07 (solid) or b20 (dotted) opacities.

The analysis by Bedin et al. (2008a) revealed a discrepancy between the cluster main sequence turn-off age and the age determined from the termination of the cooling sequence (the magnitude of the faint peak of the WD LF): the WD cluster age is lower by about 25% (the WD cluster age is about 6 Gyr, and the main sequence turn-off age is about 8 Gyr), implying a too-fast cooling of the WD models employed in the analysis (they included CO crystallization and phase separation). This discrepancy was seemingly resolved by García-Berro et al. (2010), whose WD models that also included ²²Ne diffusion had a slower cooling; they were able to recover an age of about 8 Gyr from the WD cooling sequence, which is consistent with the main sequence turn-off age.

In the intervening years, thanks mainly to the works by Brogaard et al. (2011, 2012, 2021), tight constraints have been placed on the cluster parameters from analyses of the photometry and spectroscopy of a sample of its main sequence and subgiant branch eclipsing binaries and analyses of the CMD of the cluster main sequence, red giant branch, and core He-burning phases. The resulting values for metallicity, age, reddening, and distance modulus are [Fe/H] = +0.29 ± 0.03(random)±0.08(systematic), age t = 8.3 ± 0.3 Gyr, E(B − V) = 0.14 ± 0.02, and (m − M)_V = 13.51 ± 0.06, respectively. At the same time, the very recent model calculations by Salaris et al. (2022) and Bauer (2023) show that the cooling delay induced by neon diffusion is shorter than previously calculated, leading Bauer (2023) to question whether the latest generation of WD models would still be able to lead to WD ages consistent with the cluster age determined from the main sequence.

We addressed this question by employing the recent Salaris et al. (2022) theoretical WD cooling sequences, which include the effect of neon diffusion in addition to CO crystallization and phase separation. Figure 1 displays a HST/ACS CMD of two 8.5 Gyr WD isochrones computed with Salaris et al. (2022) cooling tracks and calculated with either the Cassisi et al. (2007) electron conduction opacities (hereafter c07 opacities) or the Blouin et al. (2020) opacities (hereafter b20 opacities; see the next section for more about opacities). We used the Cummings et al. (2018) WD initial–final mass relation and progenitor lifetimes from Hidalgo et al. (2018) models¹ for [Fe/H] = +0.3, which corresponds to a metallicity Z = 0.03. We compared the WD isochrones to the cluster cooling sequence (from Bedin et al. 2008a) after being shifted by E(B − V) = 0.15 and (m − M)_V = 13.50 (consistent with the constraints highlighted above), using the transformations from A_V (for R_V = 3.1) to A_F606W and A_F814W from Bedin et al. (2005).

The F606W magnitude of the bottom end of both isochrones is clearly fainter than the observations (the isochrone calculated with the b20 opacities is the faintest one, due to the shorter cooling times of the models). This is also very clear from a comparison of the theoretical LF obtained from the isochrone calculated using the c07 opacities with the observed, completeness-corrected counterpart. This theoretical LF (as well as all the other LFs that will be discussed later in this paper) was calculated by randomly drawing progenitor masses according to a Salpeter mass function (power law with exponent equal to −2.35) and interpolating along the WD isochrone to determine the magnitude (and colour) of the corresponding WD. These magnitudes (and colours) were then perturbed by a Gaussian random error with σ values determined from the photometric analysis from Bedin et al. (2008a). The resulting number distribution of the magnitudes was then binned exactly as the observed LF.

Changing the age and distance modulus by ±0.5 Gyr and ±0.1 mag does not yet allow a match to the observed faint peak of the LF. The bottom line is that the cooling of the WD models – for both choices of opacities – seems to be too fast to match the observations, confirming the inference by Bauer (2023) that the cooling delay due to ²²Ne diffusion might not be enough to achieve consistency between WD cluster ages and the age from eclipsing binaries and the main sequence turn-off. We therefore investigated whether WD cooling tracks and isochrones that include the effect of ²²Ne distillation can lead to a better agreement between these two independent age determinations.

3. Model calculations

We computed sets of WD cooling models that include the effect of ²²Ne distillation and use the same code, physics inputs, grid of masses (M_WD = 0.54, 0.61, 0.68, 0.77, 0.87, 1.0, and 1.1 M_⊙), and initial chemical stratification (for progenitors with [Fe/H] = 0.3, which corresponds to a metallicity Z = 0.03 and hence a uniform ²²Ne mass fraction equal to 0.03 in the core) as the BaSTI-IAC WD calculations presented in Salaris et al. (2022). These calculations² include ²²Ne diffusion in the liquid phase, crystallization, and a phase separation of the CO mixture and are our baseline models for computing the cooling delay caused by the inclusion of the distillation process.

The models have a pure-H envelope that comprises a fraction q(H) = 10⁻⁴ of the total mass, M_WD, around pure-He layers of mass fraction q(He) = 10⁻², which surround the CO core, as in the baseline calculations (the standard ‘thick layers’ of H-atmosphere WD computations). As for the baseline models, we performed two sets of calculations with distillation, using the c07 and b20 electron conduction opacities, respectively. We will see that the effect of distillation on the cooling times depends quantitatively on the choice of opacities. In a nutshell, these two sets of opacities treat the regime at the transition from moderate to strong degeneracy differently; this is the relevant region covered by the H and He envelopes of the models and is still subject to sizable uncertainties (see Blouin et al. 2020; Cassisi et al. 2021; Salaris et al. 2022, for more in-depth discussions).

To include the ²²Ne distillation, we proceeded as follows. When, after a computational time step, the Coulomb parameter Γ (expressed in terms of the Coulomb parameter for a pure carbon composition, Γ_C; see below) in one or more layers reaches the critical value, Γ_cr, C, for the crystallization of the CO mixture according to the phase diagram from Blouin & Daligault (2021, see their Eq. (33) but also our Eq. (2) for a modification that accounts for the presence of ²²Ne), we checked whether the ²²Ne abundance in the surrounding liquid layers was higher than the minimum value for distillation to occur, as given by (from the calculations by Blouin et al. 2021)

$$\begin{array}{c}\hfill x_{\mathrm{Ne}}=-0.0004866{\mathrm{\Gamma }}_{\mathrm{cr},\mathrm{C}}+0.09076,\end{array}$$(1)

where x_Ne denotes the number fraction of ²²Ne nuclei.

If the current abundance was above the threshold, we started the distillation by evolving the abundances of ²²Ne, carbon, and oxygen in these layers towards the values of the distilled composition (in mass fractions) – 0.3143, 0.6857, and 0.0, respectively, which would be reached eventually when locally Γ_C = 208. The quantity Γ_C = 6^2/3e²/(a_ek_BT) denotes the value of the Coulomb parameter for a pure carbon composition, where e is the elementary charge, k_B the Boltzmann constant, a_e = (4πn_e/3)^−1/3 the mean inter-electronic distance, and n_e the electron number density. When Γ_C = 208, the distillation of neon is completed and the layers where neon has accumulated crystallize³.

While current phase diagram calculations (Caplan et al. 2020; Blouin et al. 2021) can be used to clearly identify the starting point of the distillation process (i.e. when the solid becomes lighter than the coexisting liquid) and its stopping point (i.e. when the liquid and solid phases have the same composition), the exact trajectory between these two points remains uncertain. In our calculations, we made the simple assumption of a linear variation in the composition with respect to Γ_C.

When the chemical abundances in the layers undergoing distillation are being evolved after a time step, the outer liquid layers are mixed out to the outer boundary of the core (see Blouin et al. 2021), and their abundances (uniform, because of the full mixing) were recalculated to ensure the conservation of the mass of the various elements. When distillation started at the centre, we also stopped the diffusion of ²²Ne in the liquid layers. The difference in the internal energy between the model with the new abundances and that at the previous time step is accounted for in the energy equation.

This procedure was repeated at the next time steps, and the outermost layer starting distillation moves steadily outwards with decreasing core temperature, as long as the ²²Ne abundance in the liquid phase is above the required threshold.

Figure 2 displays, as an example, the progressive change in the chemical stratification during neon distillation across the core of our 0.68 M_⊙ models calculated with the c07 opacities (the progression is the same for the models calculated with b20 opacities). After distillation starts, all chemical profiles display a flat outer part, which encompasses the liquid layers that were fully mixed during neon distillation, and an inner part, which is undergoing distillation and where the abundances display trends of increasing ²²Ne and ¹²C and decreasing ¹⁶O when moving from the inner edge of the liquid region towards the centre. The reason for this increase is that at any given time the value of Γ_C in layers progressively closer to the centre is higher (because of increasing densities and roughly constant temperatures), and hence the local abundances are closer to the values that correspond to the full distillation of neon.

Fig. 2. 12C, 16O, and 22Ne abundance profiles across 0.68 M⊙ WD models at the start of the cooling sequence (solid lines), at the end of the distillation process (dash-dotted lines), and at two intermediate stages (dotted and dashed lines in order of decreasing core temperature; see the main text for details).

As distillation progresses, the neon and carbon abundances in the liquid layers decrease and the oxygen abundance increases (as required by the conservation of the total mass of these elements). Meanwhile, the outer edge of the region undergoing distillation (which encloses a mass M_dist) moves outwards, and the inner abundances continue to evolve towards the values at full distillation.

The steadily decreasing abundance of ²²Ne in the liquid phase eventually drops below the threshold for distillation to continue (given by Eq. (1)) before Γ_C = 208 anywhere within M_dist. When this happens, distillation is stopped and M_dist reaches its final value; at the next time step of the calculations, the layers within M_dist are made to crystallize and release latent heat. Since distillation is stopped, no further compositional change is possible and therefore crystallization is inevitable. We find that for all values of M_WD, the final value of M_dist is equal to ∼25% of the M_WD.

The layers above M_dist also start to crystallize after distillation has stopped, according to the local value of the Coulomb parameter, and undergo phase separation. The CO and ²²Ne chemical profiles are now flat in these layers, with the mass fraction of ²²Ne reduced to below 0.01; the exact value of the neon abundance depends on the value of M_WD. The CO phase separation was treated according to Blouin & Daligault (2021) but with the following update to account for the presence of ²²Ne during crystallization: in this phase, we did not reactivate neon diffusion in the liquid layers, because their abundance is so low that the effect on the cooling times is negligible (see e.g. Salaris et al. 2022).

The presence of neon changes the freezing temperature of the plasma compared to the case of a binary CO mixture. To take this effect into account, we specified the crystallization coupling parameter as (Bédard et al. 2024)

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{cr},\mathrm{C}}={\mathrm{\Gamma }}_{\mathrm{cr},\mathrm{C}}^0+(c_1{x}_{\mathrm{O}}+c_2{x}_{\mathrm{O}}^2+c_3{x}_{\mathrm{O}}^3)x_{\mathrm{Ne}},\end{array}$$(2)

where ${\mathrm{\Gamma }}_{\text{cr},\text{C}}^0$ is the value given by the CO phase diagram in Blouin & Daligault (2021), x_O and x_Ne are number fractions in the liquid phase, and c₁ = 1096.69, c₂ = −3410.33, and c₃ = 2408.44. Similarly, the separation between the liquidus and the solidus, Δx_O, is modified from Eq. (34) of Blouin & Daligault (2021) as follows:

$$\begin{array}{c}\hfill \mathrm{\Delta }{x}_{\mathrm{O}}=a_0^{\prime }{x}_{\mathrm{Ne}}+(a_1+a_1^{\prime }{x}_{\mathrm{Ne}})x_{\mathrm{O}}+(a_2+a_2^{\prime }{x}_{\mathrm{Ne}})x_{\mathrm{O}}^2+a_3{x}_{\mathrm{O}}^3+a_4{x}_{\mathrm{O}}^4+a_5{x}_{\mathrm{O}}^5,\end{array}$$(3)

where $a_0^{\prime }=0.640125$, $a_1^{\prime }=2.218484$, $a_2^{\prime }=-4.599227$, and the unprimed a_i are given in Table II of Blouin & Daligault (2021). Even when the fractionation of Ne is not strong enough to lead to distillation, the solid phase is still depleted of Ne. This is taken into account using the following equation (fitted to the simulation data of Blouin et al. 2021):

$$\begin{array}{c}\hfill \mathrm{\Delta }{x}_{\mathrm{Ne}}=\min (0,-0.611587x_{\mathrm{Ne}}+0.782489x_{\mathrm{Ne}}{x}_{\mathrm{O}}),\end{array}$$(4)

where Δx_Ne is the difference in the Ne number fraction between the solid and liquid phases, and the abundances on the right-hand side of the equation are for the liquid phase. We note that Eqs. (2)–(4) were fitted to simulation data in the ranges 0.3 ≤ x_O ≤ 0.8 and 0 ≤ x_Ne ≤ 0.035 and should not be used for compositions outside this range⁴.

Figure 3 shows three snapshots of the chemical stratification for the same 0.68 M_⊙ WD models discussed before: at the start of crystallization, at an intermediate stage in the process, and when the core is fully crystallized. Crystallization does not alter the chemical abundances within M_dist, but the uniform CO and ²²Ne chemical abundances beyond the distilled core undergo a progressive change with time due to phase separation, the net effect of which is to move oxygen towards the interior and carbon towards the exterior of the core, as shown clearly in Fig. 3.

Fig. 3. Chemical abundance profiles across 0.68 M⊙ WD models at the start of crystallization (solid lines), at the end of crystallization (dashed lines), and at an intermediate stage (dotted lines). Dotted, dashed, and solid lines overlap in the inner core, which has undergone neon distillation. Dotted and dashed lines also overlap out to M/Mtot ∼ 0.6 (see the main text for details).

At the outer boundary of the distilled core region (which includes about 25% of the total WD mass), a discontinuity of the C and O abundances develops during the crystallization of the outer layers, as seen in Fig. 3. This is due to the shape of the CO phase diagram, which causes a local increase in O and a decrease in C at crystallization (see Eq. (3)) The local increase in O and decrease in C cause an inversion of the mean molecular weight, μ, and an associated maximum density discontinuity of Δρ/ρ ∼ 0.001. This value is typical for all the WD models we calculated; however, according to Blaes et al. (1990), it is not large enough to break the solid lattice in the lower-μ solid layers, which would therefore be stable to overturning.

At intermediate stages between the beginning and the end of crystallization, Fig. 3 shows that the outer layers of the core still display a uniform chemical profile (in the model shown, the crystallization front has reached layers enclosing about 60% of the WD mass). However, the abundances of neon, carbon, and oxygen are different from the values at the beginning of the process because of the abundance changes in the solid phase.

During the crystallization of the layers above the distilled core, it is unclear whether ²²Ne distillation restarts. This uncertainty is due to an oscillatory behaviour between distillation and no distillation, which in turn is due to the changes in neon and oxygen abundances in the liquid phase caused by the phase separation of the neighbouring solid layers (Eqs. (3) and (4)). Distillation would recommence when the ²²Ne abundance in the liquid layers increases above the threshold value given by Eq. (1), which is affected by the value of the abundance of oxygen in the liquid layers, and no distillation would occur when the neon abundance decreases again below the threshold. In a nutshell, phase separation tends to decrease the neon abundance in the crystallizing layers (see Eq. (4)), and as a consequence, ²²Ne increases in the surrounding liquid region, eventually becoming higher than the minimum value for distillation to start. However, when the layers just above the crystallization front start to experience neon distillation, the neon abundance in the liquid layers tends to decrease and distillation ends soon after it starts.

In some test calculations we found that if we let the distillation start again (when the conditions are met) by stopping CO crystallization and following the same procedure described before, distillation will cease almost straight away because the neon abundance in the liquid layers decreases below the values given by Eq. (1) after a few time steps. The result is that these oscillations produce narrow mass shells (thickness of the order of 1% of the total core mass) with only very small increases in the neon abundance caused by distillation and very minor contributions to the model’s energy budget.

In our final models, we assumed that distillation does not restart after it has been completed in the central layers that enclose a mass M_dist. The formation of the aforementioned narrow ‘distillating’ shells – if this happens – would in any case have no appreciable effects on our main conclusions because of their minor contributions to the energy budget.

3.1. The effect of ²²Ne distillation on cooling times

Figure 4 displays the difference in cooling times between these new calculations and the baseline models without distillation (Δt_cool) as a function of the luminosity. We show the results for sets of models calculated with the c07 opacities and b20 opacities.

Fig. 4. Difference in cooling times as a function of luminosity between WD models (masses are as labelled) calculated with and without 22Ne distillation, using either the c07 (solid lines) or b20 opacities (dotted lines).

When distillation first begins at the very centre of the models, the Δt_cool values are slightly negative because of the missing contributions from the neon diffusion in the liquid layers, and from the crystallization of the CO mixture in the layers that start distillation. This is a brief phase (it covers ∼0.1 dex in luminosity), and the Δt_cool values are never lower than ∼−0.3 Gyr.

With increasing time (decreasing luminosity), the effect of distillation overcomes the lack of energy inputs from CO crystallization and neon diffusion, and the values of Δt_cool begin to increase steeply above zero. For all masses, the Δt_cool at this point increases sharply with luminosity due to the energy gained during the distillation of ²²Ne, followed by a flatter profile after distillation has stopped. In this region, the exact trend of Δt_cool with luminosity depends on the effect of the crystallization of the CO abundance profiles above M_dist, which are different from the baseline models due to the chemical redistribution during distillation.

There is no monotonic trend in the value of Δt_cool with M_WD at the end of distillation, because the effect on the cooling times depends on the combination of how much energy is gained in the process (which increases with increasing M_WD) and the luminosity at which it is released (lower masses start distillation at lower luminosities). For a fixed amount of extra energy released at a luminosity L, the resulting Δt_coolincreases with decreasing L.

As a result, the cooling time delay caused by the distillation of ²²Ne first increases with increasing M_WD, reaching a maximum for the model with M_WD = 0.68 M_⊙. It then decreases with further increasing M_WD. Overall, the delays at the end of distillation are longer for the models calculated with the c07 opacities, which, for a given value of M_WD, start distillation at lower luminosities compared to the calculations with the b20 opacities, due to a different relationship between L and the core temperature.

4. Comparison with the NGC 6791 cooling sequence

The models described in the previous section were employed to calculate WD isochrones for both choices of opacity and various ages, using the same progenitor metallicity, lifetimes, and initial–final mass relation described in Sect. 2. Figure 5 compares the CMD of isochrones with ages of 8, 8.5, and 9 Gyr with the observed cooling sequence, after applying the same reddening and distance modulus as in Fig. 1.

Fig. 5. Same as the lower panel of Fig. 1, but showing two sets of theoretical [Fe/H] = 0.3 WD isochrones with 22Ne distillation calculated employing the c07 (upper panel) and b20 opacities (lower panel), for ages of 8.0 (dotted lines), 8.5 (solid lines), and 9.0 Gyr (dashed lines).

There are a few striking features to note in this comparison. First of all, the faint end of the 8.5 Gyr isochrones is shifted to brighter magnitudes compared to their counterparts in Fig. 1. This is the effect of the extra energy input due to the distillation of neon, which slows the evolution of the models undergoing CO crystallization.

The second obvious feature is the change in shape of the faint part of the isochrones compared to the calculations without distillation, a consequence of the variation in Δt_cool induced by distillation with M_WD. The other important feature to note is the weak effect of age (in the age range relevant to this cluster) on the CMD position and shape of the faint end of the isochrones, especially in the case of calculations with the c07 opacities. The reason is that, at the luminosities corresponding to this isochrones’ age range⁵, neon distillation is still efficient in the cores of all objects with mass between M_WD ∼ 0.6 M_⊙ and ∼0.65–0.9 M_⊙ (depending on the chosen opacities), and the extra energy gained in the process slows the cooling to such a level that their luminosity is virtually unchanged.

All these points become even clearer when we consider Fig. 6, which displays the distribution of the evolving WD mass along the isochrones of Figs. 1 and 5. The upper panel shows the increase in the brightness of all masses at the faint end of the isochrones. This is caused by the cooling delays, Δt_cool, due to neon distillation, which are listed in Table 1 for representative masses. The shift in magnitude for a given M_WD depends on the value of Δt_cool – larger for higher Δt_cool – but also on the cooling speed when distillation is not included – larger shifts for higher cooling speeds without distillation.

Fig. 6. Theoretical WD isochrones. Upper panel: distribution of the evolving WD mass as a function of the mF606W magnitude along the 8.5 Gyr isochrones with and without neon distillation. Solid black and dashed red lines denote calculations with the c07 and b20 opacities, respectively. Lower panel: same as the upper panel, but for the 8.0, 8.5, and 9.0 Gyr isochrones of Fig. 5, which all include neon distillation, calculated with either the c07 (solid lines) or b20 opacities (dashed lines).

The lower panel of Fig. 6 is similar to the upper panel, but it displays only isochrones that include neon distillation, with ages of 8.0, 8.5, and 9.0 Gyr, as in Fig. 5. It is shown very clearly that the brightness of all masses between ∼0.6 M_⊙ and ∼0.65 M_⊙ (for b20 opacities) or ∼0.9 M_⊙ (for c07 opacities) is independent of the isochrone age because these objects are still undergoing neon distillation in their cores, and their luminosity evolution essentially stalls. Higher mass objects show a variation in brightness with changing isochrone age, because they are evolving past the completion of neon distillation, and the energy contribution of this process has vanished.

This leads to the question of whether the isochrones that include neon distillation are able to reconcile the WD age with the main sequence age. Figure 7 shows, more clearly than the CMD comparison, that this is indeed the case, at least for models calculated with c07 opacities. We display here the observed LF of NGC 6791 WDs (bin size of 0.1 mag), together with the LFs calculated from our isochrones with neon distillation. We can see that in the case of calculations with the c07 opacities, the termination of the theoretical LF matches the observations. The increase in the number of objects when approaching the bottom of the theoretical LF is caused by the increasing number of WDs when M_WD decreases from 0.7–0.8 M_⊙ to ∼0.6 M_⊙ – the mass evolving at the faint end of the isochrone – due to the initial–final mass relation and the power law with a negative exponent for the chosen progenitor mass function. It is important to stress that in the case of the c07 isochrones, the magnitude of the termination of the theoretical WD LF does not change if the progenitor mass function is changed, because the m_F606W magnitude along the isochrone has an absolute maximum at M_WD ∼ 0.60 − 0.64; however, the overall shape of the LF does change.

Fig. 7. Comparison of the observed WD LF of NGC 6791 (filled circles with error bars) with theoretical LFs calculated from the same isochrones of Fig. 5. Dotted, solid, and dashed lines correspond to ages of 8.0, 8.5, and 9.0 Gyr, respectively. The normalization of the theoretical LFs is arbitrary (but the total number of stars is the same in all the theoretical LFs).

The different shape of the isochrones in Fig. 5 compared to the case without neon distillation implies that their faint end covers a very small colour range, potentially in contrast with the observed broad cooling sequence at these magnitudes. However, this small colour range is not in disagreement with the observations once photometric errors are added. This is shown by Fig. 8, which displays a synthetic CMD for the 8.5 Gyr c07 isochrone, including the appropriate photometric errors. This synthetic CMD was initially calculated as described in Sect. 2, but we then subtracted stars according to the completeness factors as a function of magnitude derived from the observations (Bedin et al. 2008a) to simulate the observed CMD displayed in Fig. 8 (the number of stars in both diagrams is the same). The bottom end of the synthetic CMD also shows a wide colour range like the observations, that is entirely due to the large 1σ error on the (F606W − F814W), which, for example, at a reference m_F606W = 28.0 amounts to ∼0.1 mag.

Fig. 8. Comparison of the observed CMD of NGC 6791 WDs with a synthetic CMD computed from the 8.5 Gyr WD isochrone (c07 opacities) shown in the top panel of Fig. 5 (see the main text for details).

The comparison with the theoretical LFs from isochrones calculated with the b20 opacities is more complex because of the different shape of the M_WD-m_F606W relation (and of the isochrones) compared to the previous c07 case. Figure 7 shows that these LFs display a too-bright faint peak compared to the observations, by about 0.2 mag. However, the situation is less clear-cut if we look again at Fig. 6. Due to the non-monotonic trend in the m_F606W with increasing WD mass in the range ∼0.6 M_⊙–∼0.7 M_⊙, the WD masses around 0.66 M_⊙ have a magnitude m_F606W ∼ 28.1 at an age of around 8.5–9.0 Gyr along the b20 isochrones; this is consistent with the magnitude of the faint peak of the observed cluster WD LF. Due to the chosen progenitor mass function (and initial–final mass relation), in the actual LF the number of objects with this mass (coming from progenitors of ∼2.5 M_⊙ according to the assumed initial–final mass relation) is lower than the number of objects with a mass of around 0.6 M_⊙ (coming from progenitors with masses equal to ∼1.5 M_⊙) that populate the peak of the LF. If the actual WD mass distribution in the cluster is instead, and somehow ad hoc, peaked around 0.66 M_⊙ – due, for example, to an initial–final mass relation, and/or a progenitor mass function different from the one assumed in these calculations, and/or dynamical evolutionary effects that have affected the mass distribution of the objects in the observed cluster field – the LF calculated with b20 opacities would show a peak consistent with the observed faint peak for ages between ∼8.5 and 9.0 Gyr.

5. An observational test with JWST

In this section we briefly discuss an observational test that could independently confirm the efficiency of neon distillation in the interiors of this cluster’s WDs, and potentially even help us identify which results – either c07 or b20 – provide a more accurate treatment of the electron degeneracy opacity.

Figure 9 displays the JWST/NIRCam (F115W, F200W) CMD of the same WD isochrones investigated so far, with and without distillation for ages of 8.0, 8.5, and 9.0 Gyr. The difference in the shape between the two sets is also remarkable in these filters. When distillation is included, the isochrones calculated with the c07 opacities display an extended (by 0.3–0.4 mag) horizontal sequence in the (F115W − F200W) colour due to the behaviour of the bolometric corrections with T_eff. This sequence is not present in the counterpart without distillation and is populated by the objects with masses between ∼0.6 M_⊙ and 0.9 M_⊙, as also shown in Fig. 9. The same objects in the isochrone without distillation are distributed along a tilted sequence in this CMD, which is also much less extended in colour than the isochrone with distillation. In the case of models calculated with the b20 opacities, the isochrones with distillation are more horizontal than those without, but less straight and with a narrower colour extension compared to the case of their c07 counterpart.

Fig. 9. WD isochrones in JWST filters. Left panels: distribution of the evolving WD mass as a function of the JWST F200W absolute magnitude along the 8.0, 8.5, and 9.0 Gyr [Fe/H] = 0.3 isochrones with neon distillation, calculated with either the c07 opacities (solid lines) or the corresponding isochrones in a JWST CMD (solid lines in the lower panel), compared to the counterparts without neon distillation (dashed lines in the lower panel). Right panel: same as the left panels, but for isochrones calculated with the b20 opacities.

Figure 10 displays the synthetic CMDs (350 stars each⁶) obtained for the 8.5 Gyr isochrones with and without distillation, shifted to the apparent magnitudes and colours expected for NGC 6791 in these filters⁷, and 1σ photometric errors in F115W and F200W equal to 0.03 mag. For a given choice of opacity, the differences between the CMDs of the two populations (with and without distillation) are still quite obvious. It is also intriguing to note that with these photometric errors, the synthetic populations calculated with distillation and the two different choices of opacity display a similar horizontal morphology at their faint, well-populated branch, but with different colour extensions. The extension in the b20 CMD is about half that of the c07 case, mirroring the differences seen between the corresponding isochrones in Fig. 9. We conclude that an observed CMD with these photometric errors should be able to further confirm (or disprove) the efficiency of neon distillation in the interiors of NGC 6791 WDs and help us identify the more appropriate treatment of the electron conduction opacity in the regime at the transition from moderate to strong degeneracy.

Fig. 10. Synthetic JWST CMDs of NGC 6791 cooling sequence. Upper panel: synthetic CMDs for 8.5 Gyr WD populations (see the main text for details) at the distance of NGC 6791, calculated from isochrones with neon distillation and either the c07 (black filled circles) or the b20 opacities (red open squares). Lower panel: same as the upper panel, but for calculations without the inclusion of neon distillation.

To assess whether these observations are feasible, we employed the most updated online JWST exposure time calculator (ETC version 3.0)⁸ and estimated a required integration time of about 5 h for each of the two filters, F115W and F200W, to attain the signal-to-noise ratio (S/N) of 25–30 needed to achieve the random photometric error used in these synthetic CMDs. For this calculation, we assumed a black-body spectral energy distribution for a temperature of 5000 K (the approximate T_eff of the faintest WDs in the simulation), normalized at m_F115W = 27.0 and m_F200W = 26.7 for the expected end of the cooling sequence. These numbers are also supported by our JWST multi-filter calibrated photometry of the open cluster NGC 2506, available from the catalogue publicly released by Nardiello et al. (2023). In those observations, a total exposure time of just 43 s in the filter F115W allowed us to reach a S/N = ∼15–30 down to m_F115W = 21. According to an approximate calculation based just on the rescaling of the exposure times, we would expect to go ∼6 mag fainter in 3 h with F115W, and similarly for the redder filter F200W.

6. Summary and conclusions

We have presented the first CO WD evolutionary models that self-consistently include the effect of neon distillation in addition to CO crystallization, phase separation and neon diffusion in the liquid phase. Our calculations included the c07 and b20 electron conduction opacities, respectively, and cover the full range of CO WD masses for a twice-solar progenitor metallicity, which is appropriate for the old open cluster NGC 6791.

We have shown that these state-of-the-art WD evolutionary models produce LFs that can match the faint end of the observed WD LF in NGC 6791, for ages consistent with the narrow range determined from eclipsing binary stars and the main sequence turn-off. Without the inclusion of distillation, the termination of the theoretical WD cooling sequences is too faint compared to the observations.

All cluster WDs with masses between ∼0.6 M_⊙ and ∼0.65 or 0.9 M_⊙ (for the b20 and c07 opacities, respectively) undergo neon distillation, and the shape of the isochrones is very different compared to the case without distillation. Also remarkable is the lack of evolution in the brightness of the LF faint end when distillation is included, which is due to the luminosity evolution of the models stalling whilst distillation is ongoing.

This means that at the metallicity and age of NGC 6791, WD age-dating on its own may not be very precise when neon distillation is included in the models’ computation, particularly for models calculated with the c07 opacities. In this case, we find that theoretical LFs with ages ranging from 7 to 10 Gyr that include the photometric errors of the cluster data all display the same magnitude of the faint peak at the termination of the cooling sequence. If we did not have tight constraints on the age of NGC 6791 from the main sequence, the WD LF – with neon distillation included – would be consistent with an age range of 4 Gyr. Even after reducing the photometric errors to just 0.02–0.03 mag across the whole LF, increasing the number of objects, and reducing the bin size by half, the magnitude of the termination of the theoretical LF is still basically the same, between 7 and 10 Gyr. This is at odds with the case of models without neon distillation, which produced LFs with varying peak magnitudes across this age range.

In the case of models calculated with the b20 opacities, this age degeneracy is greatly mitigated because only the LFs in the age range between 8.5 and 9 Gyr display the same magnitude of the peak at the termination of the cooling sequence. However, in this case, we need a fine-tuned mass distribution (peaked around M_WD ∼ 0.66) to match the observed faint peak of the cluster WD LF.

We have proposed and described new JWST observations that could independently prove the efficiency of neon distillation in the interiors of NGC 6791 WDs, based on the shape of the cluster cooling sequence in appropriately chosen JWST CMDs. An important byproduct of these proposed observations is their potential to help solve the current uncertainties on the electron conduction opacities for the modelling of WD envelopes.

Finally, as mentioned in Sect. 1, all other metals in the CO WD cores are at least one order of magnitude less abundant (in mass fraction) than ²²Ne, and their contribution to the energy budget – and hence the cooling times – might be expected to be small, if not negligible. This seems to be the case when we consider ⁵⁶Fe, the most abundant element after ²²Ne with a mass fraction equal to only about 9% of the ²²Ne abundance (see e.g. Salaris et al. 2022). Salaris et al. (2022) made use of the results from Caplan et al. (2021) regarding the phase diagram of a CO–Fe mixture to show that the sedimentation of Fe in the CO core during phase separation increases the cooling times of CO WD models by at most 100–200 Myr for a metallicity Z = 0.04 (which is even higher than that of NGC 6791), a negligible contribution compared to the effect of neon diffusion and distillation.

Our modelling work has highlighted uncertainties related to the stop and start condition for distillation that should be investigated in future work. In particular, it is unclear whether ²²Ne-rich shells can be formed after distillation is completed in the inner core. Future work should also try to improve upon the current treatment of composition change during distillation, which is currently assumed to be linear with respect to the Coulomb coupling parameter.

Acknowledgments

We thank our referee for comments and suggestions that have improved the presentation of our results. MS acknowledges support from The Science and Technology Facilities Council Consolidated Grant ST/V00087X/1. SB thanks the Banting Postdoctoral Fellowship and CITA National Fellowship programs for support and is grateful to Matt Caplan for useful discussions. SC acknowledges financial support from PRIN-MIUR-22: “CHRONOS: adjusting the clock(s) to unveil the CHRONO-chemo-dynamical Structure of the Galaxy” (PI: S. Cassisi) finanziato dall’Unione Europea – Next Generation EU, and from INFN (Iniziativa specifica TAsP).

References

All Tables

All Figures

Fig. 1. The cooling sequence of NGC 6791. Upper panel: observed WD LF of NGC 6791 (corrected for completeness) compared with a theoretical LF calculated for an age of 8.5 Gyr using a [Fe/H] = 0.3 WD isochrone (with c07 opacities) from models calculated without including the effect of 22Ne distillation (see the main text for details). Lower panel: CMD of the observed WD cooling sequence of NGC 6791 together with [Fe/H] = 0.3, 8.5 Gyr isochrones calculated without 22Ne distillation and with either the c07 (solid) or b20 (dotted) opacities.

In the text

	Fig. 2. 12C, 16O, and 22Ne abundance profiles across 0.68 M⊙ WD models at the start of the cooling sequence (solid lines), at the end of the distillation process (dash-dotted lines), and at two intermediate stages (dotted and dashed lines in order of decreasing core temperature; see the main text for details).
In the text

	Fig. 3. Chemical abundance profiles across 0.68 M⊙ WD models at the start of crystallization (solid lines), at the end of crystallization (dashed lines), and at an intermediate stage (dotted lines). Dotted, dashed, and solid lines overlap in the inner core, which has undergone neon distillation. Dotted and dashed lines also overlap out to M/Mtot ∼ 0.6 (see the main text for details).
In the text

	Fig. 4. Difference in cooling times as a function of luminosity between WD models (masses are as labelled) calculated with and without 22Ne distillation, using either the c07 (solid lines) or b20 opacities (dotted lines).
In the text

	Fig. 5. Same as the lower panel of Fig. 1, but showing two sets of theoretical [Fe/H] = 0.3 WD isochrones with 22Ne distillation calculated employing the c07 (upper panel) and b20 opacities (lower panel), for ages of 8.0 (dotted lines), 8.5 (solid lines), and 9.0 Gyr (dashed lines).
In the text

Fig. 6. Theoretical WD isochrones. Upper panel: distribution of the evolving WD mass as a function of the mF606W magnitude along the 8.5 Gyr isochrones with and without neon distillation. Solid black and dashed red lines denote calculations with the c07 and b20 opacities, respectively. Lower panel: same as the upper panel, but for the 8.0, 8.5, and 9.0 Gyr isochrones of Fig. 5, which all include neon distillation, calculated with either the c07 (solid lines) or b20 opacities (dashed lines).

In the text

	Fig. 7. Comparison of the observed WD LF of NGC 6791 (filled circles with error bars) with theoretical LFs calculated from the same isochrones of Fig. 5. Dotted, solid, and dashed lines correspond to ages of 8.0, 8.5, and 9.0 Gyr, respectively. The normalization of the theoretical LFs is arbitrary (but the total number of stars is the same in all the theoretical LFs).
In the text

	Fig. 8. Comparison of the observed CMD of NGC 6791 WDs with a synthetic CMD computed from the 8.5 Gyr WD isochrone (c07 opacities) shown in the top panel of Fig. 5 (see the main text for details).
In the text

Fig. 9. WD isochrones in JWST filters. Left panels: distribution of the evolving WD mass as a function of the JWST F200W absolute magnitude along the 8.0, 8.5, and 9.0 Gyr [Fe/H] = 0.3 isochrones with neon distillation, calculated with either the c07 opacities (solid lines) or the corresponding isochrones in a JWST CMD (solid lines in the lower panel), compared to the counterparts without neon distillation (dashed lines in the lower panel). Right panel: same as the left panels, but for isochrones calculated with the b20 opacities.

In the text

Fig. 10. Synthetic JWST CMDs of NGC 6791 cooling sequence. Upper panel: synthetic CMDs for 8.5 Gyr WD populations (see the main text for details) at the distance of NGC 6791, calculated from isochrones with neon distillation and either the c07 (black filled circles) or the b20 opacities (red open squares). Lower panel: same as the upper panel, but for calculations without the inclusion of neon distillation.

In the text