Residual Carbon in Oxygen–Neon White Dwarfs and Its Implications for Accretion-induced Collapse

Farmer et al. (2015) use MESA to investigate the properties of carbon burning in SAGB stars as a function of the initial stellar mass and rotation rate, in addition to mixing parameters, including convective overshooting. They use fine spatial and temporal resolution in order to carefully resolve the carbon flame. We therefore elect to use the MESA input files from Farmer et al. (2015) as the framework for producing ONe WDs with realistic composition profiles. We updated the Farmer et al. (2015) input files to be compatible with MESA version 9793. We then evolve models from the pre-main sequence to the SAGB and through the carbon flame phase.

We now summarize some of the key input physics and modeling assumptions made by Farmer et al. (2015) and hence adopted in our models. For complete details, we refer the reader to Farmer et al. (2015) and the MESA instrument papers (Paxton et al. 2011, 2013, 2015, 2018).

We use the MESA nuclear network sagb_NeNa_MgAl.net, which uses 22 isotopes to cover the H, He, and C burning phases (including the pp chains and the CNO, NeNa, and MgAl cycles). Mass loss is treated via the Reimers mass loss prescription (Reimers 1975) on the RGB with η = 0.5 and a Blöcker mass loss prescription (Bloecker 1995) on the AGB with η = 0.05. Convective stability is evaluated using the Ledoux criterion. Semiconvection is included via the prescription of Langer et al. (1985) with a dimensionless efficiency factor of 0.01. Thermohaline mixing is included using the prescription of Kippenhahn et al. (1980) with a dimensionless efficiency factor of 1. Convective overshooting is implemented in MESA via an exponential decay of the convective mixing diffusion coefficient beyond the fully convective boundary (Herwig 2000). The lengthscale of this decay, and hence the extent of the overshooting region, is controlled by the dimensionless parameter ${f}_{\mathrm{ov}}$, which multiplies the local pressure scale height. When we include overshooting, we use an overshooting efficiency of ${f}_{\mathrm{ov}}=0.016$. Rotation in MESA uses the shellular approximation (Meynet & Maeder 1997) and is described in Paxton et al. (2013). Stellar models are given an initial solid body rotation rate at zero-age main sequence (ZAMS). When we include rotation in our models, we use an initial rotation rate of 0.25, the critical rotation rate ${\omega }_{\mathrm{crit}}$. The MESA implementation of rotation and angular momentum transport closely follows that of Heger et al. (2000) and Heger et al. (2005). Transport is treated within the diffusion approximation and considers the effects of the following processes: the dynamical shear instability, secular shear instability, Eddington–Sweet circulation, Goldreich–Schubert–Fricke instability, and Spruit–Tayler dynamo. We do not include the effects of rotationally enhanced mass loss or other angular momentum loss mechanisms such as magnetic braking. While many of these choices have associated uncertainties and caveats, on the whole, we believe that this produces models of SAGB stars that are broadly representative of those in the literature.

Following SAGB stars through their final phases is computationally demanding. For a recent review of some of the challenges of modeling stars in this mass range, see Doherty et al. (2017). We invest the effort to carefully track the off-center carbon burning as it migrates to the center and forms the degenerate ONe core. However, difficult evolutionary phases remain before the star leaves behind a WD. One challenge comes from needing to track a large number of thermal pulses (hundreds to thousands). This would require a prohibitively large number of timesteps, so we wish to halt our models before the TP-SAGB phase.⁴ This does mean that the WD models that we produce have somewhat different masses than if they had been allowed to evolve though the TP-SAGB. Another challenge comes from contending with those models that experience "dredge-out" events. This occurs when the hydrogen-rich convective envelope merges with the helium-burning convective shell (e.g., Ritossa et al. 1999; Gil-Pons & Doherty 2010). This leads to a hydrogen flash as protons are mixed into regions of high temperature on the convective turnover time. During this phase, our MESA models encounter numerical difficulties and are unable to continue. Stellar evolution codes may be poor tools to follow these dredge-out events, as the hydrogen burning can become sufficiently rapid that standard assumptions about mixing length theory may be unlikely to hold (Jones et al. 2016a).

In summary, we cannot simply allow our MESA SAGB models to naturally evolve, eject their envelopes, and produce WDs. Since the challenges occur after completion of the carbon flame phase (and thus after the core composition is done evolving), we decide to manually halt the evolution and eject the envelope.⁵ We develop a simple, ad hoc procedure that halts the models before the onset of thermal pulses or dredge-out.

We stop the evolution and eject the envelope once two conditions are met simultaneously. First, we check if the highest temperature in the region of the star where the energy generation rate is dominated by CNO-cycle hydrogen burning is greater than roughly $3\times {10}^{7}\,{\rm{K}}$. This condition is naturally triggered as the hydrogen-burning shell develops before the first thermal pulse and is necessarily triggered before high-temperature proton ingestion can lead to numerical difficulties. Second, we check that the CO core mass and the ONe core mass are within 10% of each other. This condition ensures that we do not stop the evolution before the ONe core has finished forming. At this point, we remove the envelope using an artificially enhanced stellar wind with $\dot{M}\sim {10}^{-2}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. We turn off nuclear burning during this phase for numerical simplicity. After the envelope has been removed, we are left with a hot ONe WD model with a thin hydrogen and helium envelope.

Our goal is to produce a set of WD models that have a range of qualitatively different chemical profiles. Therefore, like Farmer et al. (2015), we vary a handful of stellar parameters including initial mass, rotation at the ZAMS, and the location and strength of convective overshooting.

We run a set of non-rotating models without any convective overshoot (series NRNO). We run a set of models (series RO) with the fiducial rotation and overshooting parameters adopted by Farmer et al. (2015), which includes overshoot at all convective boundaries. The propagation of the carbon flame, and therefore the subsequent chemical profile, is sensitive to the treatment of mixing at the convective boundary near the flame (Siess 2009; Denissenkov et al. 2013). To explore this effect, we also run a set of rotating models (series RMO) with an overshooting treatment where mixing below convective regions associated with carbon burning does not occur. This corresponds to the MESA  options

• overshoot_f_below_burn_z_shell=0.000
• overshoot_f0_below_burn_z_shell=0.000.

This choice is motivated by the findings of Lecoanet et al. (2016), who performed idealized hydrodynamic simulations of a carbon flame and found that little mixing occurs at this interface, as convective plumes were not able to penetrate into the carbon flame.

We produced a total of 12 WD models. These models and their identifiers are listed in Table 1 together with the key parameters including initial mass, rotation at ZAMS, and convective overshooting. All but one model produce massive ONe WDs. Model 5 produces a hybrid CO–ONe WD.

Table 1.  Summary of ONe WD Models

Series	Model	Initial Mass	Rotation	Overshooting	WD Model Mass	Central 12C (%)
		(M⊙)	$(\omega /{\omega }_{\mathrm{crit}})$	(${f}_{\mathrm{ov}}$)	(M⊙)	Pre-cool	Post-cool
RMO	1	7.2	0.25	0.016/0.000	1.06	6.99	2.69
	2	7.6			1.18	0.98	0.69
	3	8.0			1.28	0.22	0.18
	4	8.4			1.34	0.25	0.07
RO	5	7.2	0.25	0.016	1.08	32.72	16.30
	6	7.6			1.17	15.76	2.89
	7	8.0			1.28	3.94	1.02
	8	8.4			1.34	0.31	0.09
NRNO	9	9.5	0.00	0.000	1.01	2.46	2.65
	10	10.0			1.09	4.23	1.64
	11	10.5			1.14	2.16	1.19
	12	11.0			1.19	0.68	0.63

Note. Models in a series are evolved with an identical treatment of rotation and convective overshooting. The left part of the table shows the model identifiers. The middle part shows initial conditions and model parameters. The right part shows key properties of the formed WDs. Note that the WD progenitor models had their evolution artificially truncated before the TP-SAGB. If they were allowed to evolve though the TP-SAGB, the mass of the WD might be somewhat different (depending on the efficiency of third dredge-up).

Download table as:  ASCIITypeset image

Convective overshooting and rotation can have large effects on the chemical profile of the WD, and especially that of ¹²C. Figure 1 shows the carbon profiles in our models; each panel compares different initial masses within the same model series. Broadly, models from more massive stars have less carbon than those from lower mass stars.⁶ This trend is set as models transition from flames that fail to reach the center (leaving lots of unburned carbon) at lower mass to central carbon ignitions (with almost complete carbon burning) at higher mass. The detailed behavior in between these extremes is complex, with a mixture of carbon flames and flashes (see, e.g., Figure 3 in Farmer et al. 2015). This trend is particularly pronounced in the models that include overshooting (left and center panels). Models 1–8 from series RMO and RO have carbon depletion, especially at larger radii, while models 9–12 from series NRNO have comparatively homogeneous carbon profiles.

Zoom In Zoom Out Reset image size

Stellar models that exhibit efficient convective overshooting below the carbon flame tend to extinguish the carbon flames earlier resulting in an increased abundance of carbon in the core. This effect is clearly visible in the left and center panels of Figure 1. The profiles of series RO and series RMO, which differ only by their assumption about mixing below the carbon flame, have similar carbon profiles in the outer half of the cores. However, series RO has systematically higher carbon abundances in the inner regions than series RMO. We find that all 12 WD models exhibit qualitative differences in their carbon profiles, depending on their initial masses and modeling assumptions.

Siess (2007) constructed SAGB models using the stellar evolution code STAREVOL. These models do not take into account the effects of rotation. Like our MESA models, they use a decaying exponential treatment for convective overshooting with the same value of ${f}_{\mathrm{ov}}=0.016$, though this is only applied at the edge of the convective core. We compare our models with their solar metallicity (Z = 0.02) models that include convective overshooting.

Figure 2 shows the central ¹²C abundance at the end of carbon burning as a function of initial mass. The models from Siess (2007) are shown as gray squares. Models with higher initial masses have lower central carbon fractions. Our MESA models exhibit similar trends. Our models are not directly comparable to the Siess models at fixed initial mass because the different assumptions about rotation and overshooting imply a different mapping between initial mass and core mass. However, these independent sets of SAGB models display the same feature that stars in this transition mass range show a broad range of central carbon mass fractions ranging from tenths of a percent to tens of percent.

Zoom In Zoom Out Reset image size

Figure 2 also illustrates the sensitivity of the central carbon fraction to convective overshooting below regions associated with carbon burning. The series RO models have central carbon mass fractions up to a factor of 10 greater than those in series RMO at fixed initial mass; again, the only difference between these models is the overshooting parameter at the lower boundary of the convection zone immediately above the carbon flame.

The models in Gutiérrez et al. (2005) start with central conditions ${\rho }_{{\rm{c}}}={10}^{9}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ and ${T}_{{\rm{c}}}=2\times {10}^{8}\,{\rm{K}}$. From their Figure 3, it appears that these models almost immediately begin carbon burning. For lower carbon mass fractions, the carbon-burning runaway takes slightly longer, allowing additional compression to occur during the temperature increase. The lowest carbon mass fraction model (${X}_{12}=0.01$) appears to run out of fuel before reaching temperatures $\gtrsim {10}^{9}\,{\rm{K}}$ and then subsequently rejoins the compression-driven evolution toward higher density and the eventual ${}^{20}\mathrm{Ne}$ electron captures that are characteristic of their carbon-free models.

The fact that the models of Gutiérrez et al. (2005) appear to be burning carbon at or near their initial central conditions is puzzling. Figure 4 of Gasques et al. (2005) shows ignition lines (defined as where the energy generation rate is equal to the thermal neutrino losses) and burning timescales for pure carbon that cover the relevant range of densities and temperatures. The initial central conditions from Gutiérrez et al. (2005) have a burning timescale longer than a Hubble time.

Figure 6 shows the carbon ignition lines for a range of carbon mass fractions. These are evaluated with the same input physics as our subsequent MESA models. We use the carbon burning rate from Caughlan and Fowler (1988) and the "extended" screening treatment in MESA, as described in Paxton et al. (2011). This combines the results of Graboske et al. (1973) in the weak regime and Alastuey and Jancovici (1978) with plasma parameters from Itoh et al. (1979) in the strong regime. The neutrino loss rates are calculated via the fitting formulae of Itoh et al. (1996). For comparison, the solid gray line shows the MESA implementation of the rate given by Gasques et al. (2005) for pure carbon⁸ (${X}_{12}=1.0$). At low temperatures (roughly in the regime where T is below the ion plasma temperature), these ignition curves disagree, reflecting differences in the screening treatment. Figure 4 of Gasques et al. (2005) shows the theoretical uncertainty associated with the ignition curves in this regime; the pure carbon MESA ignition curve in Figure 6 lies within that range.

Zoom In Zoom Out Reset image size

Figure 6 also indicates where the key weak reactions will become important. The dashed lines show where the electron-capture timescale for each of ${}^{25}\mathrm{Mg}$, ${}^{23}\mathrm{Na}$, ${}^{24}\mathrm{Mg}$, and ${}^{20}\mathrm{Ne}$ is equal to a characteristic compression timescale of 10⁴ yr (Schwab et al. 2015), with more rapid electron captures occurring to the right of the line. These curves were evaluated using the rate tabulations of Suzuki et al. (2016), assuming ${Y}_{{\rm{e}}}=0.5$.

In order to illustrate the potential role of carbon, we follow the approach of Schwab et al. (2015, 2017a) by constructing quasi-homogeneous WD models with parameterized abundances. This provides a simple way to isolate the effects of ${}^{12}{\rm{C}}$. In Section 5, we will apply this understanding to the more realistic WD models described in previous sections. The calculations shown in this section use MESA version 10108. All input files will be made publicly available at http://mesastar.org.

We first create a pre–main sequence star of $1.0\,{M}_{\odot }$ with a solar composition. We select this as a starting point because a high-entropy model is more amenable to our relaxation procedure than a highly degenerate WD model. During these relaxation phases, we disable nuclear reactions. We evolve until the model reaches a central density of $\mathrm{log}({\rho }_{{\rm{c}}}/{\rm{g}}\,{\mathrm{cm}}^{-3})=4$, and then we relax its composition to one of the compositions indicated in Table 2. These are the fiducial compositions from Schwab et al. (2015, 2017a) plus a variable amount of carbon. We continue the evolution, with the model cooling and contracting, until it reaches a central density of $\mathrm{log}({\rho }_{{\rm{c}}}/{\rm{g}}\,{\mathrm{cm}}^{-3})=7$. We then have the models accrete a 50/50 oxygen–neon mixture⁹ at a rate ${10}^{-6}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ until the WD has reached $\mathrm{log}({\rho }_{{\rm{c}}}/{\rm{g}}\,{\mathrm{cm}}^{-3})\,=8.6$. This value is selected to be below the density at which electron captures on the isotopes present in the initial composition will first occur. This procedure creates near-Chandrasekhar mass WD models with homogeneous inner regions of the desired compositions.

Table 2.  The Set of Compositions Used in Our MESA Models

Isotope	SQB15C	SBQ17C
${}^{12}{\rm{C}}$	0.0–0.1	0.0–0.1
${}^{16}{\rm{O}}$	0.50	0.50
${}^{20}\mathrm{Ne}$	Remainder	Remainder
${}^{23}\mathrm{Na}$	⋯	0.05
${}^{24}\mathrm{Mg}$	0.05	0.05
${}^{25}\mathrm{Mg}$	⋯	0.01

Note. Each composition is referenced in the text by the identifier listed in the top row. Each column lists the mass fractions of the isotopes that were included. Dashes indicate that a particular isotope was not included. A range of ${}^{12}{\rm{C}}$ mass fractions was considered; the ${}^{20}\mathrm{Ne}$ mass fraction was chosen to ensure the mass fractions sum to unity.

Download table as:  ASCIITypeset image

We then take these models and continue to let them accrete at ${10}^{-6}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. This choice represents an accretion rate that is roughly characteristic of any near-Chandrasekhar WD accretor stably burning hydrogen or helium (e.g., Wolf et al. 2013; Brooks et al. 2016). The long period of accretion in these models means that the memory of their initial central temperatures are erased as the models come into a balance where the compressional heating balances neutrino losses (Paczyński 1971; Schwab et al. 2015; Brooks et al. 2016). However, this does imply that the temperature of these models is set by our choice of accretion rate. We evolve until the model experiences carbon ignition (defined as the energy release from carbon burning locally exceeding thermal neutrino losses). We note that some models, in particular those with carbon fractions ≲0.01, may evolve beyond the carbon ignition lines corresponding to their initial carbon abundances. This is because by the time an evolutionary model reaches the temperature and density of the ignition curve corresponding to its initial carbon abundance, some of the carbon can have already been consumed.

Schwab et al. (2017a) demonstrated that models that experience electron captures after significant Urca-process cooling has occurred can develop convectively unstable regions near the electron-capture front. Proper modeling of this convection likely requires both theoretical and numerical progress and remains an important open question. We discuss the difficulties of evolving MESA models in this phase more specifically in Section 5. To circumvent these issues, we artificially suppress the action of convection in the MESA models using the control mlt_option='none'. This is an important caveat that applies to the evolutionary tracks shown in Figure 7 after ${}^{24}\mathrm{Mg}$ electron captures $(\mathrm{log}(\rho /{\rm{g}}\,{\mathrm{cm}}^{-3})\gtrsim 9.6)$. Nonetheless, these models still provide insight into when carbon burning first begins to affect the evolution.

Zoom In Zoom Out Reset image size

Figure 7 shows the evolution of the central temperature and density in the models. The thin lines show the results of models with the SQB15C composition, which do not include Urca-process cooling, as there is no ${}^{23}\mathrm{Na}$ or ${}^{25}\mathrm{Mg}$. These models are less physically realistic, but we show them primarily for comparison with the results of Gutiérrez et al. (2005). Models with ${X}_{12}\gtrsim 0.02$ experience carbon ignition at densities below where electron captures on ${}^{24}\mathrm{Mg}$ would occur. Energetically, this requires burning approximately $0.02\,{M}_{\odot }$ of carbon to heat the WD to the conditions for dynamical burning (Piro & Bildsten 2008). Therefore, because these models have roughly that much carbon, if we continued their evolution, they could reach the conditions for thermonuclear explosions. Note that the ${X}_{12}=0.1$ model is similar to what one would expect from a fully mixed hybrid CO–ONe WD. Models with very low carbon abundances ${X}_{12}\lesssim 0.003$ evolve much like carbon-free models; the thin orange and blue lines are nearly identical and continue until oxygen ignition. Intermediate abundance models with $0.003\lesssim {X}_{12}\lesssim 0.02$ experience carbon ignition as a result of the A = 24 electron captures. This range encompasses the carbon fractions in many of the ONe WD models shown in Section 3.

However, the presence of Urca-process cooling changes this picture, causing models across a range of compositions and accretion rates to reach ${T}_{{\rm{c}}}\lesssim {10}^{8}\,{\rm{K}}$ (see Equation (15) in Schwab et al. 2017a). The thick lines in Figure 7 show the results of models with composition SBQ17C; the presence of ${}^{23}\mathrm{Na}$ and ${}^{25}\mathrm{Mg}$ leads to significant cooling. Models that experience this cooling will therefore almost certainly reach the ${}^{24}\mathrm{Mg}$ electron captures before igniting carbon. Models with very low carbon abundances ${X}_{12}\lesssim 0.003$ evolve much like carbon-free models (though again, we note that these models do develop convectively unstable regions that are artificially not allowed to mix). Models with ${X}_{12}\gtrsim 0.01$ experience carbon ignition as a direct result of the A = 24 electron captures.

Stellar models indicate significant abundances of the odd mass number isotopes (in particular, ${}^{23}\mathrm{Na}$) that are responsible for the Urca-process cooling. Thus the SBQ17C models are the ones that should guide our understanding of WDs with self-consistently set compositions. Urca-process cooling prevents carbon ignitions from occurring in advance of the exothermic electron captures on ${}^{24}\mathrm{Mg}$. The SBQ17C models do show carbon ignition triggered by the A = 24 electron captures above a critical carbon fraction of around 1% by mass. Since this is within the range seen in our self-consistent stellar models of ONe WDs, this is an indicator that carbon may play an important role in at least some systems.

In the overall evolutionary scenario for AIC, the ONe WD will stop cooling once it begins accreting from its companion via Roche lobe overflow. This can be achieved through expansion of the secondary driven by stellar evolution processes or a shrinking of the orbit through angular momentum losses. Either process results in a massive WD accreting material.

We model the effects of accretion from a companion by adding a constant accretion rate of ${}^{16}{\rm{O}}$. We again select a rate of ${10}^{-6}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ as roughly characteristic of any near-Chandrasekhar WD accretor stably burning hydrogen or helium. Accreting oxygen simplifies the models, as it avoids having to follow the numerically challenging hydrogen, helium, or carbon flashes that could occur during the accretion process. We turn rotation off for all models to avoid numerical problems with angular momentum transport in the WD. This is consistent with past work that has primarily studied non-rotating WDs.

We use a version of the original nuclear network sagb_NeNa_MgAl.net modified to add the isotopes ${}^{19}{\rm{O}}$, ${}^{20}{\rm{O}}$, ${}^{20}{\rm{F}}$, ${}^{21}{\rm{F}}$, ${}^{23}{\rm{F}}$, ${}^{23}\mathrm{Ne}$, ${}^{24}\mathrm{Ne}$, ${}^{24}\mathrm{Na}$, ${}^{25}\mathrm{Ne}$, ${}^{27}\mathrm{Na}$, and ${}^{27}\mathrm{Mg}$, and the weak reactions that link them to the isotopes originally present in our network. We use the reaction rate tables from Suzuki et al. (2016), which are specifically calculated for the conditions in high-density ONe cores. We adopt EOS options¹⁰ and spatial/temporal resolution controls similar to those of Schwab et al. (2017a), as it is demonstrated that those choices lead to numerically converged models (see their Appendix B).

We show the evolution of the central density and temperature of a subset of our models in Figure 8. In particular, we select models 2 and 6 from Table 1. These models both descend from 7.6 ${M}_{\odot }$ ZAMS models, with the only difference between them being the efficiency of overshooting associated with the carbon flame. These produce similar WD models that differ mainly in their carbon abundance; model 2 and model 6 have central carbon mass fractions of 0.007 and 0.03, respectively.

Zoom In Zoom Out Reset image size

Compared to the toy models in Section 4, these models begin somewhat cooler. This is because the toy models had been steadily accreting since they were much lower mass and had reached a balance between compressional heating and neutrino cooling (Paczyński 1971). In contrast, our realistic ONe WD models had already cooled below this temperature and, being initially more massive, have not accreted for long enough to reheat.

The annotations on the lower part of the figure indicate the location of key weak reactions. We see some difference from the toy models shown in Figure 7. Mildly exothermic captures on ${}^{27}\mathrm{Al}$ occur; this species is not included in our toy models. The Urca-process cooling due to ${}^{25}\mathrm{Mg}$ is less than in the toy models, due both to the initially lower temperature and the lower abundance of this isotope. The cooling from ${}^{23}\mathrm{Na}$ is more dramatic, though the small glitch near the end of cooling is an artifact of the reaction rate tables and does correspond to slightly too little cooling (see Appendix D in Schwab et al. 2017a). The separation (in density) of the ${}^{24}\mathrm{Mg}$ and ${}^{24}\mathrm{Na}$ electron captures is caused by different assumed strengths of the (currently unmeasured) non-unique second forbidden transitions (see Section 5 in Schwab et al. 2017a).

All model tracks are similar up until the electron captures on ${}^{24}\mathrm{Na}$. The heating from these is sufficient to ignite carbon. We are unable to evolve the models, including the effects of convection, beyond this point. By artificially suppressing convection, we are able to continue their evolution. These tracks are marked "no convection" in Figure 8. Subsequently, the difference between the two model tracks is apparent, with the a bifurcation reflecting the difference in the carbon abundance. The low carbon fraction model 2 continues to high density, while the low carbon fraction model 6 experiences a relatively low-density carbon-assisted oxygen ignition.

Both in Section 4 and in this section, we encountered numerical difficulties when convective regions developed in the regions of the star undergoing electron-capture reactions. We circumvented this problem by artificially turning off convection. The further evolution of the models is thus unphysical, though hopefully remains schematically useful.

To understand why it is difficult to evolve the models beyond this point, it is first important to recall why these models become convective. The onset of exothermic electron captures leads to a strongly destabilizing temperature gradient. However, this is locally accompanied by the development of a strongly stabilizing composition (electron fraction) gradient. When using the Ledoux criterion, this can lead to overall convective stability (e.g., Miyaji & Nomoto 1987; Schwab et al. 2015). The crossed temperature and composition gradients mean that this region is unstable to double-diffusive convection, though this process is thought to operate too slowly to have a significant effect, given the short evolutionary timescale. When including the effects of Urca-process cooling, the models of Schwab et al. (2017a) became convectively unstable because thermal conduction causes a more rapid spread of the thermal gradient and leads to a situation in which the thermal and compositional gradients can no longer balance locally because their spatial extents are no longer well-matched. We emphasize that the development of these convective regions is unrelated to the presence of ${}^{12}{\rm{C}}$ and occurs even in carbon-free models.

However, once the exothermic electron captures raise the temperature enough for carbon burning to begin, its energy release can further contribute to convective instability.¹¹ The development of a convective core has important consequences in the overall evolution of the star. Most importantly, it distributes the entropy generated by continued carbon burning and exothermic electron captures over a larger mass. This has been shown to lead to oxygen ignition at higher densities (relative to convectively stable models), favoring collapse to an NS (Miyaji et al. 1980; Miyaji & Nomoto 1987).

The numerical difficulties that occur once a convection zone develops are essentially the same as those encountered in previous work; thus, we briefly recapitulate the discussion from Section 6 of Schwab et al. (2017a).

When using MLT in MESA, cells that are convective have roughly the adiabatic temperature gradient, while non-convective cells have the radiative temperature gradient. Therefore, when cells change between convective and radiative, they change the overall temperature profile, and in turn, this can alter the convective stability of neighboring cells. (This is further complicated by the fact that the sub-threshold electron-capture reaction rates are exponentially sensitive to temperature.) As MESA iterates to find the solution to its equations, convective boundaries change location iteration-to-iteration and the solver fails to converge.

It is clear that regions become convectively unstable, but an understanding of how and whether these convective regions grow remains an outstanding problem. If a long-lived central convection zone does form, matters are further complicated, as this allows for the operation of the convective Urca process (Paczyński 1973). This continues to be difficult to model in stellar evolution codes (e.g., Lesaffre et al. 2005), and current versions of MESA produce physically inconsistent results during such phases (Schwab et al. 2017b).

The choice to artificially suppress convection corresponds to the assumption that these convective regions do not grow. The evolution of our "no convection" models will be roughly accurate only if the the energy and composition changes from nuclear reactions remain approximately local. In this case, the presence of carbon leads to ignition at a factor of $\approx 2$ lower density than in the carbon-free case. In the case that a large convective core does form (invalidating our local treatment), it seems clear that this leads to higher density oxygen ignition than in a case without a convective core; in the carbon-free models of Gutierrez et al. (1996), non-convective models ignite at $\mathrm{log}({\rho }_{{\rm{c}}}/{\rm{g}}\,{\mathrm{cm}}^{-3})\approx 10$, with ignition in convective cores occurring at a central density a factor of $\approx 2$ higher. However, the presence of carbon is likely to reduce the density relative to a carbon-free model, as the additional heating source means that fewer electron captures are required. Since we cannot follow the models through a convective Urca phase, we are unable to reliably estimate the size of this shift.