Laminar Flame Speeds in Degenerate Oxygen-Neon Mixtures

The collapse of degenerate oxygen-neon cores (i.e., electron-capture supernovae or accretion-induced collapse) proceeds through a phase in which a deflagration wave ("flame") forms at or near the center and propagates through the star. In models, the assumed speed of this flame influences whether this process leads to an explosion or to the formation of a neutron star. We calculate the laminar flame speeds in degenerate oxygen-neon mixtures with compositions motivated by detailed stellar evolution models. These mixtures include trace amounts of carbon and have a lower electron fraction than those considered in previous work. We find that trace carbon has little effect on the flame speeds, but that material with electron fraction $Y_e \approx 0.48-0.49$ has laminar flame speeds that are $\approx 2$ times faster than those at $Y_e = 0.5$. We provide tabulated flame speeds and a corresponding fitting function so that the impact of this difference can be assessed via full star hydrodynamical simulations of the collapse process.


INTRODUCTION
Degenerate oxygen-neon (ONe) cores with masses near the Chandrasekhar mass can form in the evolution of ≈ 8 − 10 M single stars (e.g., Miyaji et al. 1980;Miyaji & Nomoto 1987), in interacting binary systems with varying degrees of envelope stripping (e.g., Tauris et al. 2015;Poelarends et al. 2017), in binary systems with an accreting ONe white dwarf (WD) (e.g., Canal & Schatzman 1976;Nomoto & Kondo 1991), or as the result of the merger of two WDs (e.g., Saio & Nomoto 1985;Brooks et al. 2017). When the core reaches a central density ≈ 10 10 g cm −3 , exothermic electron captures on 20 Ne occur and lead to the initiation of a deflagration wave ("flame") that propagates outward. The competition between the energy release from this flame and the electron-capture reactions on its ashes determines whether this leads to an explosion (resulting in partial or total disruption of the star) or implosion (resulting in the formation of a neutron star (NS)). This situation has long been known to be finely balanced (e.g., Nomoto & Kondo 1991;Canal et al. 1992), though the general conclusion by the end of the 1990s was in favor of collapse to a NS. Recent multidimensional simulations have reiterated that the outcome is sensitive to modeling choices and reopened the possibility that at least some cases may lead to a thermonuclear explosion (possibly also leaving a low-mass bound remnant) instead of collapse to an NS (Jones et al. 2016;Leung et al. 2019;Jones et al. 2019). One of the key ingredients in this modeling is the speed at which the flame propagates. Timmes & Woosley (1992), hereafter TW92, calculated the physical properties of conductively-propagated laminar burning fronts in high-density, degenerate carbon-oxygen (CO) and ONe mixtures. We repeat a similar set of calculations using Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al. 2011Paxton et al. , 2013Paxton et al. , 2015Paxton et al. , 2018Paxton et al. , 2019, but extend these results to a wider range of compositions motivated by expectations from detailed models of the internal composition of ONe WDs (e.g., Garcia-Berro et al. 1997;Iben et al. 1997;Siess 2006). Section 2 describes the methods we use to calculate the laminar flame speeds. Section 3 reproduces the TW92 results in both CO and ONe mixtures. We then focus on the laminar flame speeds in ONe mixtures under different conditions. In Section 4 we show how the flame speeds are only mildly affected by the presence of small amounts of 12 C, but in Section 5 demonstrate the significant influence of the electron fraction (Y e ) of the material. Section 6 provides tabulated flame speeds and a corresponding fitting function. Section 7 briefly describes the implications for models of the collapse of ONe cores.

METHODS
We use MESA revision r12115 (Paxton 2019). The input files necessary to reproduce our work are publicly available 1 and an illustration of this capability has been included in the test suite case conductive flame.
We create an initial, spatially-uniform MESA model with a temperature T b = 3 × 10 8 K, specified density ρ 9 = ρ/(10 9 g cm −3 ), and a specified unburned composition. These properties characterize the cold material into which the flame will propagate. So long as the upstream temperature is much less than the downstream (post-burn) temperature of ≈ 10 10 K, the temperature jump across the flame is approximately the same and the initial temperature does not play an important role. The total mass M defines the size of the (Lagrangian) computational domain. Because the flame width λ varies substantially with density, our domain size must as well. We always ensure that M/(ρλ 3 ) 1, but typically choose this ratio to be ∼ 100 to limit the computational cost. This also implies that λ r, so the flame is effectively planar. In practice, M ∼ 10 − 10 5 g. The small size of the domain implies that the pressure gradient due to gravity is negligible. The inner boundary is at r = 0. The outer boundary has a fixed temperature T b and a fixed pressure equal to the initial pressure of the material.
We then insert a hot spot at the center with fractional size in mass q s = 3 × 10 −4 at a temperature T s = 8 × 10 9 K (for CO mixtures) or T s = 10 × 10 9 K (for ONe mixtures). This hot spot should have a size of order the flame width and a temperature of order the postburn temperature to ensure a steady-state flame. If the hot spot is too small the flame will die. If the hot spot is too large the flame will exhibit a de-accelerating transient. With a well chosen spot size much smaller than the domain, once the flame has propagated a few flame lengths, the initial condition will be effectively erased.
We define the location of the flame to be the location of the maximal rate of nuclear energy release in the  Figure 1. Schematic of a steady-state flame. The fuel in the upstream material is initially heated by conduction until the temperature becomes large enough to ignite nuclear reactions. At the critical temperature, the energy generation rate equals the energy conduction term (heating equals cooling). The downstream material burns to its nuclear statistical equilibrium state. The entire structure, approximately isobaric, propagates into the upstream fuel with a unique speed and width.
domain (i.e., peak nuclear heating, see Figure 1). We evolve the model until the flame has propagated through 90% of the domain and then extract the steady-state flame properties. By repeating this process for different initial conditions, we calculate the laminar flame speed as a function of ρ 9 and composition. In Appendix A, we demonstrate that our results are insensitive to the details of the initial conditions and are numerically converged.

Microphysics
As discussed in TW92, the flame will have a width such that the diffusion timescale across it is comparable to the timescale at which nuclear reactions heat the material. This argument leads to an estimate of the flame speed, where D th is the thermal diffusion coefficient, nuc is a characteristic specific rate of energy generation from nuclear reactions, and e is a characteristic specific energy. Therefore, the speed of the flame is set by the energy generation rate as determined from the nuclear network and the assumed thermal transport properties of the degenerate plasma. In MESA the thermal conductivity of the degenerate electrons, accounted for via a conductive opacity (κ cond ∝ D −1 th ), comes from tables privately communicated by A.Y. Potekhin (see Section A.3 in Paxton et al. 2013). TW92 describe in detail their implementation of the transport properties; a substantially similar approach is adopted in Timmes (2000), hereafter T00. The source code for the transport properties assumed in these works is publicly available. 2 Figure 2 compares the conductive opacities over the range of temperatures and densities considered in this paper. The left panel shows the T00 and MESA values of κ cond and the right panel their relative difference. Both sources show similar density scalings, but the T00 values scale less steeply with temperature such that while the values agree at 10 9 K, κ cond is higher by a factor ≈ 2 at 10 8 K and lower by a factor ≈ 1.5 at 10 10 K. Since these variations are not systematically in the same direction, their effect is difficult to estimate, but given the scaling in Equation (1), variations in κ cond at this level correspond to ≈ 30% variations in the flame speed.

Nuclear reaction rates
The currently applicable default inputs for nuclear reaction rates are described in Appendix A.2 of Paxton et al. (2019). Rates are taken from a combination of NACRE (Angulo et al. 1999) and the Joint Institute for Nuclear Astrophysics REACLIB library (default version, dated 2017-10-20) (Cyburt et al. 2010). The MESA screening corrections are from Chugunov et al. (2007), which includes a physical parameterization for the intermediate screening regime and reduces to the familiar weak Graboske et al. 1973) and strong (Alastuey & Jancovici 1978;Itoh et al. 1979) limits at small and large values of the plasma coupling parameter.
Relatively large nuclear networks are required to fully capture the energy generation rate in neutron-rich compositions at these thermodynamic conditions. TW92 illustrate the increase in flame speed with increasing network size (their Table 5) and adopt a 130 isotope network. We perform a similar exercise, using 3 built-in MESA networks (approx21, mesa 204, mesa 495), a network constructed with the same elements as the TW92 130 isotope network (see their Table 1), and also an adaptive network that automatically adds and removes isotopes and which settles in at around 320 isotopes. Figure 3 shows the flame location as a function of time for a set of runs for a fiducial CO mixture (panel a) and a fiducial ONe mixture (panel b). Networks of more than 200 isotopes appear to be required before network size no longer makes an appreciable difference in the flame speed. This result is consistent with Chamulak et al. (2007), hereafter C07, who found that for flames in CO mixtures a 430 isotope network gave speeds up to ≈ 25% greater than a 130 isotope network. We run with 495 isotopes unless otherwise stated.
The JINA REACLIB polynomial fits to the reaction rate data end at 10 10 K as do the tabulated partition functions used to calculate the reverse rates and ensure detailed balance. Above 10 10 K, MESA fixes the rates to be their T = 10 10 K values. In some cases, especially for the ONe flames, the temperature exceeds 10 10 K and the peak in nuc occurs near the temperature thresh- 21 isos (140) 130 isos (160) 204 isos (191) 495 isos (197) adaptive (197) (a) CO mixture with X C = 0.5 at ρ 9 = 6.  old. If the MESA treatment underestimates the true peak of nuc , then this can lead to an underestimate of the flame speed. (For example, if we truncate the rates at T = 8 × 10 9 K, the ONe flame in Figure 3 has a speed of 82 km s −1 , a ≈ 10% reduction.) However, we have physical reasons to expect that this effect is not particularly large. By ≈ 1.2 × 10 10 K photodisintegration is strong enough to decompose nuclei into neutrons, protons, and alpha particles. This is an endothermic process, meaning there is an upper limit to how much more positive nuc can be achieved beyond the place where MESA truncates the rates.
Our results depend slightly on our adopted rate sources. If we use pure JINA REACLIB defaults (eliminating NACRE), the flame speeds increase. For the calculation shown in Figure 3, the result with the 495 isotope net and only the JINA rates is 207 km s −1 for the CO case and 94.1 km s −1 for the ONe case. These represent an approximately 5% speed up.
Thus there is some systematic uncertainty from nuclear reaction rates in our results which is difficult to characterize, but seems unlikely to be smaller than ≈ 10%. We note that both the above caveats result in even faster flame speeds than the ones we will report. The open and reproducible nature of our work allows this problem to be easily revisited, enabling the impact of future experimental and theoretical progress in the relevant reaction rates to be quickly assessed.

COMPARISON WITH PAST WORK
First, we consider CO mixtures. Following TW92 we select a 12 C mass fraction and put the remainder in 16 O. Figure 4 compares our results with those of TW92. Qualitatively, the agreement is good, and we reproduce the trends with ρ 9 and X C . Quantitatively, above ρ 9 = 4, our results are ≈ 5 − 10% slower, while below ρ 9 = 4, our results are faster, up to ≈ 40% at ρ 9 = 1. Figure 4 also compares the subset of our results that overlap with C07. We agree well at X C = 0.5. (It is difficult to see the symbols as they overlap.) We note, as do C07, that their fitting function does not appear to do a good job of matching their tabulated results. Our results are slower at X C = 1, though we note the C07 values also disagree with TW92 and that the primary focus of C07 was on X C = 0.3 − 0.7.
Next, we consider ONe mixtures. Following TW92, we first select a mass fraction X O of 16 O. When the mixture is not pure 16 O, we also choose a mass fraction 0.1 of 24 Mg. The remainder is 20 Ne. Figure 5 compares our results with those of TW92. There is qualitative agreement, with a trend (as in the CO case) that our flame speeds are faster than TW92 below ρ 9 = 4, up to ≈ 50% at ρ 9 = 1. Above ρ 9 = 4 and for X O = 0.6 and X O = 1.0, the agreement is within ≈ 10% of TW92. For X O = 0.8, the agreement is somewhat worse and the speeds are systematically ≈ 15% lower above ρ 9 = 6. This section demonstrates that our speeds are generally in good agreement (≈ 10%) with the results of past work. Relative to TW92, our calculations adopt a larger nuclear network (leading to faster flames), but have slightly higher conductive opacities (leading to slower flames). In the end, these effects may offset some-2 4 6 8 10 This Work Figure 4. Comparison with TW92 for CO mixtures. We compare our results (solid circles) to both their tabulated speeds (lighter, outlined symbols) and provided fitting function (line). We do the same for the subset of conditions that have data from C07 (Chamulak et al. 2007). This is tabulated points and a fit for XC = 0.5, ρ9 ≤ 6 and points for XC = 1.0, ρ9 ≤ 4. what. We have no reason to expect exact agreement with TW92 as this is not an identical calculation.

INFLUENCE OF TRACE CARBON
ONe cores are formed after off-center carbon ignition occurs and a convectively-bounded carbon deflagration propagates to the center (e.g., Farmer et al. 2015). Incomplete carbon burning that occurs as the flame approaches the center can leave residual carbon mass frac- tions of up to a few percent. Schwab & Rocha (2019) performed calculations of accreting ONe WDs including the presence of this carbon and concluded that models are unlikely to reach carbon ignition (and subsequently oxygen ignition and the formation of the deflagration) below the threshold density for 24 Mg electron captures. Here we explore whether, once the deflagration is ignited, the carbon affects the flame speed.
We select the X O = 0.6, X Ne = 0.3, X Mg = 0.1, composition used by TW92 and add a small amount of 12 C, reducing the 16 O mass fraction accordingly. Figure 6 shows the ratio of this flame speed to the carbon-free speed shown in Figure 5. The flame speed increases, reflecting the additional energy release from fusion of 12 C (relative to the 16 O that it replaced). However, for carbon mass fractions of a few percent, the flame speed increases only by ≈ 10%. Therefore, we conclude that the presence of small amounts of carbon is unlikely to have a significant effect on the laminar flame speeds.

INFLUENCE OF LOWER ELECTRON FRACTION
Detailed models of ONe WDs do not give compositions that are only 16 O, 20 Ne, and 24 Mg. Several neutronrich isotopes are typically present at mass fractions of a percent or more, meaning that material is expected to have Y e significantly below the Y e = 0.5 value of a 16 O/ 20 Ne/ 24 Mg mixture. Table 1 summarizes the abundances in the ONe core of a representative stellar model from Siess (2006). This mixture has Y e ≈ 0.49. See their Section 5 for an explanation of this core nucleosynthesis.
As the core slowly grows and its density increases further, the Fermi energy of the degenerate electrons rises. Electron-capture reactions on a given isotope become energetically favored when material exceeds its threshold density. 3 In Table 1, we indicate isotopes where these electron captures are likely to occur before the formation of the oxygen deflagration (meaning that their a This model has a central density ≈ 7 × 10 7 g cm −3 . Indicated isotopes can undergo additional electron captures as the density increases towards ≈ 10 10 g cm −3 , the density at which the oxygen deflagration is expected to form. threshold densities are below the threshold density of 20 Ne, which is ≈ 10 10 g cm −3 ). In what follows, we focus on the most abundant of these isotopes, 23 Na and 24 Mg. The effective threshold density for 23 Na is ≈ 1.6 × 10 9 g cm −3 and for 24 Mg is ≈ 4 × 10 9 g cm −3 . The timescales for the electron capture reactions are typically shorter that the evolutionary timescale of the object, so they are expected to completely convert the parent isotope to its daughter. The electron captures imply that Y e spatially varies through the core, with Y e becoming lower at higher density. By time the deflagration forms and begins to propagate, electron captures have already completely converted the 23 Na to 23 Ne and the 24 Mg to 24 Ne over the inner ≈ 0.2 M of the star. 4 For the mixture in Table 1, this further reduces Y e to ≈ 0.485.
In a concluding comment, TW92 note that Y e < 0.5 is expected and mention two calculations including reduced Y e in the form of 22 Ne. They report that for a flame in CO with Y e ≈ 0.498 the speed decreased by 4 A representative Ye profile as a function of mass is shown in Figure 11 of Schwab et al. (2017). One caveat is that if a large core convection zone were to develop, as happens in models adopting the Schwarzschild criterion for convection (e.g., Miyaji et al. 1980), the central region would likely be homogenized. ≈ 5% and for a flame in ONe with Y e = 0.48 the speed decreased by ≈ 30%. For CO flames, the effect of 22 Ne was studied by C07. They found the opposite sign of the effect, with a 22 Ne mass fraction of 0.06 leading to a ≈ 30% increase in the flame speed.
To quantify the effect of lower Y e , we calculated flame speeds at ρ 9 = 10 with a variable amount of neutron-rich material. We performed a set of calculations using each of 22 Ne, 23 Ne, and 24 Ne. In all cases, the mixture had a mass fraction 0.6 of 16 O with the remaining material being 20 Ne. Figure 7 shows the significant impact of the neutron richness, with the flame speed relative to that at Y e = 0.5 doubling by Y e ≈ 0.488. The sequences with the different isotopes overlap, indicating the speedup is largely independent of the neutron source.
The small change in Y e does not significantly affect the internal energy or thermal conductivity, but does lead to a significant change in nuc as the initial source of extra neutrons opens additional energy producing reaction channels. We examined the peak values of nuc in the calculations shown in Figure 7 and confirmed that the increasing flame speed is due to an increasing nuc at lower Y e and that it quantitatively follows the expectation from Equation (1).
To illustrate that this implies a density-dependent enhancement of the flame speed over the TW92 result, we construct two sets of models that initially have a mass fraction X 23 = 0.05 of A = 23 elements and X 24 = 0.05 of A = 24 elements. In one, the A = 23 material is always 23 Na and the A = 24 material is always 24 Mg. In the other, the spatially-uniform composition is selected differently depending on the chosen ρ 9 . The A = 23 material is 23 Na if the density is below its threshold density and 23 Ne if it is above it, while the A = 24 material is 24 Mg if the density is below its threshold density and 24 Ne if it is above it. We then run these models and extract their flame speeds. Figure 8 compares these two sets of calculations by showing the ratio of the flame speed in the case where the initial material has electron captured to the case where it has not. Above both threshold densities, where Y e has fallen to ≈ 0.49, the flame is ≈ 80% faster.

FITTING FORMULA
To allow this important effect to be incorporated in hydrodynamics calculations, we provide a simple fitting function like that of TW92, but including Y e as an additional parameter. As shown in Figure 7 Figure 9 plots the results.
At Y e = 0.5, the flame speed always increases with increasing oxygen abundance (as found in TW92). However, in our results at lower Y e , this is no longer true. Incorporating the effect of X O in the fit would require something beyond the power-law scaling used in the fit of TW92. Given the relatively weak dependence on X O , we circumvent this complication and propose the following simple fitting function that includes only ρ 9 and Y e : As shown in Figure 10, the fit agrees with the calculated points within 10% relative error at ρ 9 > 4, with the maximum error growing to a 30% underestimate at ρ 9 = 1. This fitting function will also do a worse job in pure oxygen mixtures (a relative error ≈ 30% for the X O = 1 points shown in Figure 5), but such pure mixtures are unlikely to arise in astrophysical contexts. If a more precise reproduction of our results is desired, the flame speed values are provided in Table 2, allowing direct interpolation in our results.

SUMMARY AND CONCLUSIONS
Using MESA calculations that resolve the structure of conductively-propagating deflagrations, we calculated laminar flame speeds in oxygen-neon mixtures. These speeds are a necessary ingredient in simulations of the final stages of electron-capture supernovae and accretioninduced collapse.
We demonstrated that the values of Y e ≈ 0.48 − 0.49 expected in these objects lead to an increase in the flame speed by a factor of ≈ 2 over that at Y e = 0.5, the value assumed in the widely-used prescription of Timmes & Woosley (1992). The low Y e is due to the nucleosynthesis during the helium and carbon burning phases that preceded the formation of the ONe core and to subsequent electron captures on isotopes initially present in the ONe core (most importantly 23 Na and 24 Mg) that occur as the core grows. As shown in Figure 8, this implies that the realized enhancement is density-dependent and most significant for ρ 4 × 10 9 g cm −3 (i.e., above the 24 Mg threshold density).
Full star hydrodynamics simulations that follow the propagation of the deflagration through the ONe core do so by including a sub-grid model for the flame. These models enhance the laminar speed by including a subgrid model of the flame-turbulence interaction (which allows for a larger, non-planar area to undergo burning), such that the laminar speed is only a lower limit. Eventually, this speed becomes irrelevant once the turbulence is fully developed, as a turbulent deflagration no longer depends on the laminar speed. In Section 6, we provide a tabulated set of laminar flame speeds as well as a convenient fitting function. These are suitable for incorporation into sub-grid flame models. Jones et al. (2016) show the laminar and turbulent flame velocities in their 3D hydrodynamic simulations. Typically, these flames remain laminar for ≈ 0.4 s, corresponding to ≈ 100 km of flame propagation. Typically, the inner ≈ 200 km is above the 24 Mg threshold density and thus at the lowest Y e . Therefore, we believe the factor of 2 speedup is representative of what will be realized in practice. The more rapid release of energy associated with a faster flame pushes models in the direction of being more likely to explode (meaning less likely to form a neutron star). The full implications of our results await the incorporation of this updated prescription in multidimensional models.

Figures 4 and 5 in
We thank Sam Jones, Pablo Marchant, Bill Paxton, and Fritz Röpke for helpful conversations. We thank the referee for a useful report. This work benefited from the May 2019 Lorentz Center program Electron-Capture-Initiated Stellar Collapse. We are grateful to Zoë Weber-Porter for performing some exploratory calculations of this problem as part of her UC Santa Cruz undergraduate thesis. We acknowledge use of the lux supercomputer at UC Santa Cruz, funded by NSF MRI grant AST 1828315. We thank Brant Robertson for the rapid and friendly technical support that enabled us to make use of this machine. This research was supported by the National Science Foundation (NSF) under the Software Infrastructure for Sustained Innovation program grants (ACI-1663684, ACI-1663688, ACI-1663696 Software: MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018 In this appendix, we demonstrate that the flame speeds we report are only weakly dependent on the details of the initial conditions and the spatial and temporal resolution of the MESA calculations. Figure 11 shows this for a flame in a CO mixture and Figure 12 does so for a flame in an ONe mixture. The discussion below applies equally to both figures. Table 2. Flame speed (in km s −1 ) for models described in Section 6 and shown in Figure 9.

XO
Ye ρ9 Panel (a) illustrates the effect of varying our procedure for extracting the flame speed. By default, we measure the flame speed using the change in position over the second half of the time interval covered by the calculation (indicated as [0.5, 1.0] in the legend). So long as we avoid the transients during the early part (roughly the first quarter) of the calculation, the extracted velocities are consistent at the percent level.
Panel (b) and panel (c) illustrate the effect of varying the temperature of the initial hot spot T s and its fractional size q s . As expected, so long as the hot spot causes a steady-state flame to propagate, these choices have no effect on the flame speed.
Panel (d) illustrates the effect of varying the upstream temperature T b . The flame speed increases with increasing upstream temperature, but such that a factor of ≈ 2 change in T b leads to only 1 − 2% increase in the flame speed. We would expect this to remain true so long as T b is much less than the post-burn temperature of ≈ 10 10 K.
Panel (e) illustrates the effect of varying the spatial resolution of the MESA calculation. MESA adaptively refines its mesh based on a set of mesh functions. The maximum cell-to-cell variation of these functions is maintained at around the value of the control mesh delta coeff which is set equal to 1 in our calculations. One of the built-in mesh functions has the form T function1 weight × log(T /K). This function ensures that temperature gradients are resolved, placing approximately T function1 weight zones per dex change in temperature. The number of zones in the calculation (which is ≈ 1000 for the default) varies roughly linearly with T function1 weight. The results are approximately independent of the spatial resolution, with a sub-percent increase between the default and higher resolution cases.
Panel (f) illustrates the effect of varying the temporal resolution of the MESA calculation. The control varcontrol target limits the fractional step-to-step variation of quantities in the same cell. The number of timesteps in the calculation (which is ≈ 2000 for the default) varies roughly linearly with the inverse of varcontrol target. The results are approximately independent of the time resolution, with a roughly 1% increase between the default and the highest resolution.