Neon Cluster Formation and Phase Separation during White Dwarf Cooling

The Gaia space observatory has determined parallax distances to large numbers of Galactic stars (Babusiaux et al. 2018), which allow for unprecedented tests of white dwarf (WD) models and evolution. Massive WDs have high central densities and strong gravities. Furthermore, some of them may have formed via mergers (Hollands et al. 2020) and they may have interesting compositions. Recently, Cheng et al. (2019) found that the population of massive WD known as the "Q branch" appear to have an additional heat source that maintains a luminosity of order 10⁻³ L_⊙ for gigayears. Latent heat from crystallization (Tremblay et al. 2019; Winget et al. 2009; Horowitz et al. 2010) and gravitational energy released from conventional ²²Ne sedimentation (Bildsten & Hall 2001; Deloye & Bildsten 2002; García-Berro et al. 2008; Hughto et al. 2010) do not appear to be large enough to explain this luminosity (Camisassa et al. 2020; Cheng et al. 2019). Heating from conventional electron capture and pycnonuclear (or density-driven) fusion (Salpeter & van Horn 1969; Yakovlev et al. 2006; Horowitz et al. 2008) reactions appear to need even higher densities and may depend too strongly on the density and or temperature (Horowitz 2020). There are many works discussing dark matter interactions in WDs (see, for example, Bertone & Fairbairn 2008; Graham et al. 2018, 2015; Hurst et al. 2015; Bramante 2015; Acevedo & Bramante 2019). However, a WD may be too small to capture enough dark matter for its annihilation to produce the necessary heat (Horowitz 2020).

The gravitational potential of a WD is large. Therefore, there is possibly enough energy available from sedimentation of neutron-rich ²²Ne to provide the necessary heating (Bildsten & Hall 2001). However, the expected diffusion constant for ²²Ne in a C/O mixture (Hughto et al. 2010) is too small to allow enough sedimentation before the C/O mixture freezes. In addition, ²²Ne sedimentation is significantly slowed down by C/O crystallization. Therefore, even though sedimentation is a large enough energy source, in practice, sedimentation is likely slow and most of this energy may remain untapped by the time the star freezes (Deloye & Bildsten 2002; García-Berro et al. 2008; Hughto et al. 2012, 2011).

It is possible that the ²²Ne sedimentation rate is enhanced because of the formation of mesoscopic ²²Ne clusters. Recently, Blouin et al. (2020) speculated that Ne phase separation produces mesoscopic clusters that could enhance the heating from conventional Ne sedimentation. Numerical calculations by Bauer et al. (2020) consider this possibility in more detail, showing that clustering of N nuclei enhances the downward drift rate proportional to N^2/3. This suggests that the appropriate diffusion timescale is achievable with long-lived clusters of only few hundred to a few thousand ²²Ne ions, and is validated with MESA models.

In this Letter we directly address this possibility with molecular dynamics (MD) simulations of Ne microcrystals in C/O mixtures and new calculations of the C/O/Ne phase diagram. We find that our Ne microcrystals are unstable in C/O liquid mixtures. Therefore, enhanced heating from cluster formation or phase separation is unlikely to be important in conventional C/O/Ne mixtures.

IUMD is a CUDA-Fortran classical molecular dynamics code and has been extensively used to model astromaterials in WDs and neutron stars (Caplan & Horowitz 2017). In our model nuclei are fully ionized and treated as point particles of charge Z_i that interact via a screened two-body potential

$$\begin{eqnarray}&&V({r}_{{ij}})=\displaystyle \frac{{Z}_{i}{Z}_{j}{e}^{2}}{{r}_{{ij}}}\exp (-{r}_{{ij}}/\lambda ).\end{eqnarray} \tag{ 1 }$$

with periodic separation r_ij and screening length ${\lambda }^{-1}\,=2{\alpha }^{1/2}{k}_{F}/{\pi }^{1/2}$ using ${k}_{F}={(3{\pi }^{2}{n}_{e})}^{1/3}$. The screening assumes relativistic electrons, which is not the case in WDs. At 10⁶ g cm⁻³, k_F ≈ 0.4 MeV, so we slightly overestimate λ; the true λ is 21% smaller. At 10⁷ g cm⁻³, the true λ is 7% smaller. This is expected to have little impact on our MD (see Hughto et al. 2011; Hamaguchi et al. 1997). We use λ/a = 2.816 for our mixture (consistent with Hughto et al. 2010), which forms a body-centered cubic (bcc) lattice at low temperature (Vaulina et al. 2002).

The mixture we consider has number abundances $\vec{x}\,=({x}_{{\rm{C}}},{x}_{{\rm{O}}},{x}_{\mathrm{Ne}})\approx (0.49,0.49,0.02)$ and is well motivated astrophysically. The ⁴He(2α, γ)¹²C and ¹²C(α, γ)¹⁶O reactions set x_C/ x_O ≈ 1, with order 10% variation due to reaction rates and WD mass (Lauffer et al. 2018). During helium burning the remaining CNO elements are thought to burn via ${}^{14}{\rm{N}}{(\alpha ,\gamma )}^{18}{\rm{F}}{({\beta }^{+})}^{18}{\rm{O}}{(\alpha ,\gamma )}^{22}\mathrm{Ne};$ assuming solar metallicity, x_Ne ≈ 0.02 in the WD (Bildsten & Hall 2001).

The crystallization of a one-component plasma (OCP) depends on the dimensionless parameter ${{\rm{\Gamma }}}_{i}={Z}_{i}^{5/3}{e}^{2}/{a}_{e}{k}_{B}T$ with ion charge Z_i, electron separation ${a}_{e}=(3/4\pi {n}_{e})$, and thermal energy k_BT. The OCP is solid (liquid) above (below) Γ_crit ≈ 175 (Potekhin & Chabrier 2000). Screening may raise ${{\rm{\Gamma }}}_{\mathrm{crit}}$ to about 178 (Vaulina et al. 2002). Averaging over mixture components for the multi-component plasma (MCP) gives ${{\rm{\Gamma }}}_{\mathrm{MCP}}=\langle Z{\rangle }^{5/3}{e}^{2}/{a}_{e}{k}_{B}T;$ crystallization is more complicated and sensitive to the exact mixture. A binary C/O mixture has Γ_MCP ≈ Γ_crit ≈ 230, while C/O/Ne mixtures with x_Ne = 0.2 have been found with MD to have 275 ≲ Γ_crit ≲ 300 (Hughto et al. 2011). If Ne clusters are unstable then we expect that the ${{\rm{\Gamma }}}_{\mathrm{crit}}$ we find will be closer to 300 than 230, as the region of the phase diagram probed will be that of high Ne concentration at the C/O liquid-Ne cluster surface.

Our MD mixture contains 8000 ¹²C, 8000 ¹²O, and 384 ²²Ne ions (N = 16384 = 2¹⁴, for GPU threading). Our initial configuration is shown in Figure 1 and consists of a neon microcrystal embedded in a C/O liquid in a cubic volume with periodic boundaries. The neon crystal is prepared by trimming ions from the edges and corners of a cubic bcc crystal to produce a truncated cuboctahedron containing 253 ions, 6a on its longest diameter. The C/O fluid and the remaining ²²Ne have random initial positions around the microcrystal; if the microcrystal is stable we expect it to grow by adsorption of the remaining 131 ²²Ne.

Zoom In Zoom Out Reset image size

We run isothermal simulations to resolve ${{\rm{\Gamma }}}_{\mathrm{crit}}$ and study the stability of the neon cluster at a range of temperatures. Constant temperature is approximately achieved by rescaling the velocities to a Maxwell–Boltzmann distribution with the desired temperature every 100 time steps. Our simulations therefore do not conserve energy; instead, in the long time limit our simulations trend toward equilibrium so E(t) allows us to resolve melting or freezing (e.g., heats of fusion). Our simulations differ from past work (e.g., Hughto et al. 2011) and our asymptotic states may not be true equilibrium (i.e.,  they may be superheated/cooled) because we only run as long as needed to verify stability or instability of the neon cluster.

In Table 1 we list isothermal simulations run using the initial configuration shown in Figure 1. These simulations were run for between 10⁵ and 10⁶ MD time steps with dt = 1/18ω_p with ion plasma frequency ${\omega }_{p}={(4\pi {e}^{2}\langle Z{\rangle }^{2}n/\langle M\rangle )}^{1/2}$. We clearly resolve a first-order transition between 287 < Γ < 289. At Γ ≥ 289 the C/O liquid is supercooled and quickly begins crystallizing, first nucleating around the Ne cluster (Figure 2, top left) before growing to fill the volume (top right). At Γ ≤ 287 we find that the Ne microcrystal dissolves into the liquid (bottom left) and in the long time limit is fully mixed in the volume (bottom right). At Γ = 244 the microcrystal immediately melts to a liquid and mixes, while at higher Γ the cluster seems to sublimate as the core remains solid while ions slowly desorb into the liquid. When the configuration with solid C/O is evolved for longer times we also resolve slow diffusion on the lattice, which seems to disperse the Ne-enriched core. We have thus shown that our initial condition is unstable, suggesting that ²²Ne microcrystals are not expected in C/O WDs.

Zoom In Zoom Out Reset image size

Our Γ_crit ≈ 288 is high compared to past work for C/O/Ne mixtures, which find Γ_crit ≈ 230 (see Table 1 in Hughto et al. 2012). Our system may be finding the melting point of a very Ne-rich system. Furthermore, the simulation does not have enough time for diffusion to bring the composition of the solid phase into equilibrium with the composition of the liquid phase. Both finite-size and finite-time effects may be important in comparing to an equilibrium phase diagram (computed below).

We now consider what impurities might form stable microcrystals. Generally speaking, stronger separation is observed in mixtures with greater contrast in charge Z. While ²³Na or an isotope of Mg may be present in comparable abundances to ²²Ne they may not have large enough Z to strongly separate. Simulations with Mg in place of the Ne at Γ = 290 also show a slow sublimation of the cluster into the liquid, though on slightly longer timescales than the Ne. In simulations with Fe run at Γ = 240, 262, and 289 the microcrystal persists after a few times 10⁷ time steps. A configuration evolved for 4.1 × 10⁷ time steps at Γ = 289 is shown in Figure 3. In all three simulations the microcrystal shows exchange with the Fe in the background and evolution in morphology. Ions on raised facets seem more likely to escape or migrate to adjacent faces to produce larger smooth surfaces, possibly an octahedron (stability, growth, and diffusion of various microcrystal morphologies may be of interest to future authors). Thus, while ²²Ne does not form stable microcrystals, higher-Z impurities are viable candidates for phase separation and clustering.

Zoom In Zoom Out Reset image size

Although mixtures with x_Ne ≈ 0.02 cannot form stable microcrystals, mixtures with larger Ne abundances could. Therefore, we compute the ternary phase diagram to determine what Ne abundance may be required for such strong phase separation. We use the code⁴ developed in Caplan et al. (2018) that implements the semianalytic method of Medin & Cumming (2010). This method identifies pairs of points on the minimum free energy surfaces that share a tangent plane (i.e.,  the double tangent construction); these points correspond to coexisting solid and liquid compositions. Compositions lying on the tangent line connecting them are therefore unstable and phase separate. For a detailed discussion see Medin & Cumming (2010) and Caplan et al. (2018).

In Figure 4 we show C/O/Ne phase diagrams for three Γ. We report temperature in Γ_C ∝ 1/T. Orange liquidus points are connected to corresponding blue solidus by green tie lines corresponding to the tangent in free energy. Below (above) the orange (blue) curve is stable liquid (solid), while the green tie lines span the unstable region.

Zoom In Zoom Out Reset image size

At the lowest temperature (Figure 4(a)), it is possible to form solid particles that are substantially enriched in Ne. The phase diagram shows solid–solid coexistence with strongly Ne-enriched mixtures, where both the Ne-enriched and Ne-depleted crystals show x_Ne > 0.3 (for readability we exclude solid–solid tie lines). We see the emergence of stable liquid near the C/Ne axis (bottom), found in the white region under the liquidus near x_C ≈ 0.75. Along the C/Ne axis we see solid–liquid phase coexistence with strong separation; given Z_Ne/Z_C = 1.66 we resolve eutectic separation.

At intermediate temperature (Figure 4(b)) we reach ${{\rm{\Gamma }}}_{\mathrm{crit}}$ for our mixture, $\vec{x}$ = (0.49, 0.49, 0.02), which is found near the middle of the right axis. In the unstable region we find separation largely consistent with the known C/O phase diagram from MD (see Figure 5 in Hughto et al. 2012). Separation is approximately independent of the Ne fraction for x_Ne ≲ 0.3. Observe that the coexistence lines are approximately parallel to the C/O-axis, which implies they fall on lines of constant Ne. Therefore, the Ne fraction in the solid and liquid are nearly equal, which suggests that any C/O/Ne mixture with x_Ne ≲ 0.3 is incapable of separating to form any Ne-rich components that could then preferentially settle out. This is qualitatively consistent with past MD (see Hughto et al. 2011) that studied mixtures of 0.02 ≤ x_Ne ≤ 0.20.

We also resolve a small set of mixtures that undergo three-component eutectic separation at intermediate Γ, similar to the mixtures with large Z₂/Z₁ and Z₃/Z₁ studied by Caplan et al. (2018). The central "wedge" between the break in tie lines contains unstable mixtures that do not fall on a single coexistence line so they cannot separate into a single solid and single liquid composition, but they can separate by forming appropriate amounts of the two solids and one liquid at the corners of this region.

At high temperature (Figure 4(c)) the solidus and liquidus curves are continuous and the region for the eutectic separation has closed. Though mixtures with x_Ne ≈ 0.30 now show weak enhancement in solidus Ne, those with x_Ne ≲ 0.30 are still largely consistent with two-component C/O separation without any Ne enrichment or depletion.

Even if there is factor of 2 variation in x_O/x_C it only translates our mixture parallel to the O-axis, which does not move our mixture into a region where it achieves strong Ne purification at any Γ. We conclude that the behavior of ²²Ne microcrystals in MD is fully consistent with the known C/O/Ne phase diagram, and that no mesoscopic effects exist that make ²²Ne cluster formation likely.

As past MD has only studied phase coexistence up to 20% Ne, we also perform a few MD simulations to validate the behavior around 30% Ne. Similar to Hughto et al. (2012) and Caplan et al. (2018), configurations are prepared by joining a cubic bcc crystal on one face with a cubic volume of liquid (as in Figure 1 in Hughto et al. 2012 and Figure 1 in Caplan et al. 2018). Compositions for the solid and liquid are chosen to approximately match predictions from the phase diagram at Γ_C = 183, allowing us to verify the strong break in the center of Figure 4(b). Our first mixture has low Ne, $\vec{x}$ = (0.4, 0, 3, 0.2), which separates to a solid ${\vec{x}}_{s}$ ≈ (0.4, 0.3, 0.3) and a liquid ${\vec{x}}_{l}$ ≈ (0.6, 0.2, 0.2). Our second mixture has high Ne, $\vec{x}$ = (0.32, 0, 15, 0.53), and strongly separates to ${\vec{x}}_{s}$ ≈ (0.1, 0.15, 0.75) and ${\vec{x}}_{l}$ ≈ (0.55, 0.15, 0.3). These simulations were evolved for 10⁷ time steps with little observed evolution. Runs varying Γ_C up and down respectively find quenching of diffusion in the liquid and melting of the crystal, suggesting these temperature variations have moved our mixtures into the region of stable liquid and stable solid. Taken together, these simulations have qualitative agreement with our phase diagram, and we conclude that the separation behavior near 30% Ne concentration is likely physical. Future work may be interested in performing a more thorough survey at high Ne concentration with MD, though this may have limited astrophysical relevance.

We find that ²²Ne microcrystals are always unstable in a C/O liquid. Either the temperature is high enough that the crystal melts and the Ne dissolves into the liquid, or the whole system including the C/O mixture freezes. Note that even at temperatures below the melting point of pure Ne, but above the C/O melting point, a large entropy of mixing causes the small concentration of Ne to dissolve into the bulk liquid. The C/O/Ne phase diagram suggests that much more Ne is necessary before it phase separates. One needs not 2% but ≈30% or more. As a result, a conventional C/O WD with n_Ne = 0.02 is not expected to form stable neon clusters with enhanced sedimentation. In summary, we find that there are no conditions where a ²²Ne-enriched cluster is stable in a C/O WD, and therefore enhanced diffusion of ²²Ne cannot explain the Q branch.

What compositions could then explain the heating that Cheng et al. infer? As seen in our phase diagram, unless the C/O ratio or ²²Ne abundance is tuned to extremes we do not expect strong ²²Ne separation, so we suggest that Q branch WDs may have an anomalous composition. For example, Camisassa et al. (2020) have suggested that x_Ne = 0.06 can provide heating on the desired timescale considering only single-particle diffusive settling rather than clusters. Another possibility is ≈1% abundance of another impurity, besides ²²Ne, with an even larger charge Z that would allow it to phase separate even when Ne does not. This impurity would need to be neutron-rich (Z/A < 0.5) to be a sedimentary heat source and have an abundance of 1% or more for there to be enough gravitational energy available. Our MD with a high-purity microcrystal shows that Z = 11 Na and Z = 12 Mg should not strongly separate in a C/O mixture. Phase diagrams of C/O/Na and C/O/Mg mixtures produced using our semianalytic method (omitted for length) are similar to the C/O/Ne in that they do not separate to form a solid enriched in the high-Z impurity when it is only abundant at the percent level, so isotopes such as ²³Na or ²⁶Mg are poor candidates for clustering.

Iron-group elements provide another possibility, as we readily observe long-lived Fe microcrystals in MD. As Γ_Fe ≈ 10Γ_C/O, Fe in C/O will phase separate and likely does not require fine-tuning of the mixture. While sedimentation of 0.1% Fe by mass may produce some notable heating, if some astrophysical process enriches Q branch WDs up to ≈1% mass fraction then settling out of Fe could provide heating for several Gyr as an otherwise conventional C/O WD cools. Thus, this work motivates including Fe in WD cooling models. This will require new phase diagrams of Fe and a survey with MD of the clustering and the characteristic sizes of Fe clusters, which will be the subject of future work.

We thank S. Cheng and E. Bauer for helpful discussions. C.J.H.'s research was supported in part by US Department of Energy Office of Science grants DE-FG02-87ER40365 and DE-SC0018083. The authors acknowledge the Indiana University Pervasive Technology Institute for providing supercomputing and database, storage resources that have contributed to the research results reported within this Letter. This research was supported in part by Lilly Endowment, Inc., through its support for the Indiana University Pervasive Technology Institute. A.C. is supported by an NSERC Discovery Grant, and is a member of the Centre de recherche en astrophysique du Québec (CRAQ).