Dimming the Lights: 2D Simulations of Deflagrations of Hybrid C/O/Ne White Dwarfs Using FLASH

The dimmest and most numerous outlier of the Type Ia supernova population, the Type Iax event, is increasingly being found in the results of observational campaigns. There is currently no single accepted model to describe these events. This 2D study explores the viability of modeling Type Iax events as a hybrid C/O/Ne white dwarf progenitor undergoing a deflagration using the multiphysics software FLASH. This hybrid was created using the stellar evolution code MESA, and its C-depleted core and mixed structure have demonstrated lower yields than traditional C/O progenitors in previous deflagration-to-detonation studies. To generate a sample, 30 “realizations” of this simulation were performed, the only difference being the shape of the initial match head used to start the deflagration. Consistent with earlier work, these realizations produce the familiar hot dense bound remnant surrounded by sparse ejecta. Our results indicate that the majority of the star remains unburned (∼70%) and bound (>90%). Our realizations produce total ejecta yields on the order of 10−2–10−1 M ☉, ejected 56Ni yields on the order of 10−4–10−2 M ☉, and ejecta kinetic energies on the order of 1048–1049 erg. Compared to yields inferred from recent observations of the dimmest Type Iax events—SN 2007qd, SN 2008ha, SN 2010ae, SN 2019gsc, SN 2019muj, SN 2020kyg, and SN 2021fcg—our simulation produces comparable 56Ni yields but too-small total yields and kinetic energies. Reignition of the remnant is also seen in some realizations.


INTRODUCTION
Type Ia supernovae (SNe Ia) are "standardizeable" candles -their lightcurves align after applying the Phillips relation (Phillips 1993), which relates their maximum brightness and decline rate.This hallmark makes SNe Ia a tape measure for the universe, a property that has been used to demonstrate that the universe is expanding at an accelerating rate (Riess et al. 1998;Perlmutter et al. 1999) and to provide a more precise calculation of cosmological parameters to test dark energy theories (Abbott et al. 2019), among others.The accuracy of these distance measurements, and therefore the accuracy of their applications, depends on how well understood the explosion mechanism is; not all SNe Ia are created equal.
SNe Ia are thermonuclear explosions of a white dwarf (WD), a burned-out star that is subject to a mass limit, the Chandrasekhar mass (M Ch ).Upon obtaining extra mass from a companion star, the WD reignites fusion reactions and subsequently explodes if it is composed principally of C/O.The mass limit of a WD greatly limits the amount of material that can be ejected from the system, including the radioactive 56 Ni whose decay gives the event its brightness.
Although the resulting light and elemental composition of a supernova are able to be observed, both the progenitor model and the explosion mechanism remain elusive.These are the keys to understanding their diversity, which indicates that the initial WD (progenitor) does not always have the same mass or elemental composition, and may have different processes causing its explosion (explosion mechanism).Especially interesting are the outliers of this population, which do not follow the Phillips relation.The dimmest and most numerous of these outliers are categorized into their own class: Type Iax events.
This paper is organized as follows.Section 1 describes the background and characteristics of SNe Ia and their dim subclass Type Iax (Section 1.1), as well as the proposed explosion mechanism (Section 1.2) and progenitor system (Section 1.3) used for this study.A brief review of related previous work is given in Section 1.4.We then describe our methodology in detail in Section 2, including the creation of our progenitor model (Section 2.1); simulation of the explosion (Section 2.2 and Appendix A); and formation of a sample (Section 2.3).Section 3 shows how we calculate the yields of our simulations.The results are shown in Section 4 and their implications are discussed in Section 5. Finally, Section 6 describes the next steps we will take to perform our future 3D study.

Type Iax Events
A subset of SNe Ia were first classified under the name Type Iax by Foley et al. (2013).Also called "02cx-like" events as a homage to their characteristic member, this group encompasses 5 -31% of all SNe Ia, and contains over 50 members, of which the most recently analyzed is SN 2021fcg (Karambelkar et al. 2021).Type Iax events can be distinguished from standard SNe Ia by their lower and wider range of peak magnitudes (-19 < M V,peak < -13), reduced ejecta velocities (2000 < v ej < 8000 km/s) and explosion energies (∼ 10 49 − 10 51 ergs), slower decline rates in the redder band lightcurves, the absence of a second maximum in the near infrared band lightcurve, significant ejecta mixing, slow spectral evolution, and divergence from the Phillips relation (Jha 2017).Type Iax events also produce much smaller ejecta yields, inferred from observations to be on average 0.5 M ⊙ , ∼ 0.003 -0.3 M ⊙ of which is 56 Ni (Foley et al. 2013) (see Section 4.8 for observational yields).As these characteristics make Type Iax events substantially different than SNe Ia, it has been suggested that their explosion mechanisms are different.As described in the next section, some leading models for a SN Ia involve a combination of a deflagration (subsonic) and a detonation (supersonic), while for Type Iax events, the hypothesized explosion mechanism is a deflagration alone.

Explosion Mechanism: Deflagration
Deflagrations have been suggested as an explosion mechanism that could reproduce the dimness and small ejected masses of Type Iax (See Section 1.1) events (Foley et al. 2013;Jha 2017).A deflagration evolves by conduction-induced burning and is therefore subsonic, meaning that the flame propagates slowly enough that the star has time to react, typically by expansion.This expansion reduces the fuel density, which causes the burning to extinguish before consuming the entire star.The result is a slower, less violent process than either a detonation alone or a deflagration-todetonation transition (DDT), both of which involve supersonic burning that can consume the star in a few tenths of a second.The gentler deflagration mechanism is expected to produce much less iron-group element (IGE) dominated material, and therefore much less 56 Ni, than a detonation or a DDT, and will accordingly produce a dimmer event.In similar studies, this deflagration explosion mechanism has been referred to as a 'failed' (Foley et al. 2013) or 'pure' (Leung & Nomoto 2020;Lach et al. 2022) deflagration.
One approach to creating a deflagration in a WD is the single-degenerate picture, in which a WD accretes mass from a companion star until thermonuclear runaway occurs.The composition and mass of the companion as well as the accretion method are variable, but the combination must produce an accretion rate that allows the WD to grow in mass, either by creating conditions for stable burning or by means of multiple nova explosions that leave some of the accreted mass behind.Interestingly, the one and only progenitor system that has ever been observed, completely serendipitously, for a Type Iax event is that of a He-donor companion star for SN 2012Z, giving credence to this hypothesis (McCully et al. 2014(McCully et al. , 2022)).The single-degenerate picture has not been able to explain the full variation in Type Iax events.An alternative being explored is the double-degenerate picture, where a deflagration is ignited by the accretion disk of a WD-WD merger (Karambelkar et al. 2021).There is also observational evidence in support of the double-degenerate picture: the only observed Type Iax event remnant, the infrared source 2MASS J00531123+6730023 or "Parker's star", has been argued to be the result of a C/O + O/Ne WD-WD merger (Gvaramadze et al. 2019;Oskinova et al. 2020;Ritter et al. 2021).This study's focus is placed on the singledegenerate picture.
The end state of simulations of deflagrations often shows a two-part system; as the deflagration is not strong enough to unbind the entire star, a large part of the star remains bound and is mixed with burning ashes, while the rest of the mass is energetic enough to be ejected from the system (Kromer et al. 2015;Leung et al. 2015;Lach et al. 2022).Section 1.4 gives an overview of the most recent of these studies in the single-degenerate, centrally ignited, near-Chandrasekhar mass scenario.Simulations of traditional C/O WD progenitors undergoing deflagrations produce an amount of 56 Ni consistent with high and intermediate brightness Type Iax events.One suggested method to reach the dimmest events is to instead use a hybrid C/O/Ne WD progenitor.

Hybrid Progenitors
White dwarfs are inactive cores -the end state of low mass (≲ 8 -10 M ⊙ ) stars.Their lighter outer layers have either dissipated as a planetary nebula or, if the star is in a binary pair, been pulled away by their companion.What is left behind is a dense core supported by electron degeneracy pressure -this is a white dwarf (WD).The WD's elemental composition is governed by the mass of the initial star, which determines how far down the chain fusion can occur.Stars with masses < 8 M ⊙ form carbonoxygen (C/O) WDs, while heavier ones (8 -10.5 M ⊙ ) create oxygen-neon (O/Ne) WDs.These C/O WDs will produce a SNe Ia under the right accretion conditions, while O/Ne WDs will collapse before ignition can occur.(Hansen et al. 2004).
Recent research has postulated a new "hybrid" C/O/Ne model -a WD made of primarily 12 C, 16 O, and 20 Ne whose interior has homogenized due to mixing.These hybrids are thought to form when the WD's carbon burning is quenched by convective boundary mixing before it reaches the core , creating a C-rich C/O core surrounded by a O/Ne shell (Denissenkov et al. 2013).It has been demonstrated that these layers mix as the WD cools due to the electron-to-baryon ratio (Y e ) gradient between the heavier O/Ne/Na mantle and the lighter C/O core (Brooks et al. 2017).The burned material has a smaller Y e due to the beta de-cay of 13 N that occurs during the burning.This gradient becomes unstable to convection as the WD cools: as heat is taken away from the mantle through neutrino transport and conduction, the mantle starts to sink and mix with the core material.This convective mixing reduces the Y e gradient, which moves to a uniform value as the system becomes isothermal.Thermohaline convective mixing continues the process until the entire interior of the star is fully mixed and the Y e gradient is 0, producing a flat composition profile.
The relatively short cooling times allow for the WD's interior to be fully mixed before the onset of accretion.This mixing significantly decreases the carbon fraction in the core.This homogenization due to mixing has been demonstrated in both 1D stellar evolution calculations by Brooks et al. (2017) and 2D direct numerical simulation by Schwab & Garaud (2019).
In this paper, we refer to a "hybrid WD" as a C/O/Ne WD that has completely mixed during cooling.The hybrid progenitor used for this study is described in detail in Section 2.1 and has been shown to produce lower 56 Ni and ejecta yields in DDT simulations due to its lower carbon fraction (Willcox et al. 2016;Augustine et al. 2019).

Previous Deflagration Studies
There are several studies that perform full hydrodynamic simulations of deflagrations of WDs in the centrally-ignited near-Chandrasekhar mass scenario.The availability of computational resources that satisfy the intense requirements of such simulations, as well as the complexity of developing software that reproduces the behavior of a thermonuclear flame on an unresolved scale, are two of the main roadblocks in achieving this goal.Over the past decade, great strides have been made in both areas, and it is now possible to perform systematic studies that investigate not only nucleosynthetic yields and observables, but also the effects of the initial flame ignition.LEAFS (Reinecke et al. 1999a(Reinecke et al. ,b, 2002) ) and FLASH (Fryxell et al. 2000;Calder et al. 2002) are two software systems capable of performing such 3D simulations, and a third by Leung et al. (2015) is currently functional up to 2D studies.
Here, we outline results from the most recent studies of deflagrations in the literature from the three codes: Lach et al. ( 2022 2022) is the newest study to date, and gives an excellent contemporary overview of the state of the field.We briefly describe the current state of research here.

C/O WDs
The more traditional C/O WD progenitor has been simulated undergoing a deflagration by multiple studies.Here, we present the most recent studies that use each of the three codes.Long et al. (2014) performed deflagrations of C/O WDs using an earlier version (Calder et al. 2007;Townsley et al. 2007) of the software used for this study, FLASH (see Section 2.2).Their 1.365 M ⊙ progenitor was made of 50/50 C/O, and ignited by placing varying numbers (63 -3500) of 16 km burned spherical regions within a central sphere of 128 -384 km radius.The end state of the simulation is recognized as having a dense core surrounded by less dense material, however unlike most other studies, all of this material is considered to be ejecta.It is likely, especially since they observe material sinking towards the core, that only some of this material is actually ejected, while the majority recombines to form a bound remnant.This is corroborated by the fact that the generated light curves do not match those of observations, though the simulation with the fewest number of ignition points is the closest.The amount of 56 Ni in the core vs. the surrounding material is related to the number of ignition points used; simulations with fewer ignition points show the majority of 56 Ni in the less dense surrounding material, while those with more ignition points show the majority of 56 Ni trapped in the core.Lach et al. (2022) perform 3D simulations of deflagrations of C/O progenitors using LEAFS.It utilizes a parabolic piecewise method (PPM) hydrodynamics solver, a FFT-based gravity solver that solves the full Poisson equation, a subgrid turbulence model, and Timmes & Arnett (1999)'s Helmholtz EOS with Coulomb corrections.The flame front is tracked using a level-set method, and nuclear burning is approximated by a reduced reaction network of 5 pseudospecies (α, 12 C, 16 O, IME, and IGE), where the density of the cell determines the burned composition.The NSE state is dynamic and is adjusted due to electron capture.It uses two nested expanding grids to simultaneously track the remnant and the ejecta.Elemental abundances were calculated from the temperature and density history of Lagrangian tracer particles added to the simulation using the nuclear network code YANN, and light curves were generated using the radiative transfer (RT) code ARTIS.Lach et al. (2022) use a traditional equal parts C/O progenitor with a central density of 2.6 × 10 9 g/cm 3 , central temperature of 6 × 10 8 K, and solar metallicity.They ignite it in a single spot, using 8 overlapping spheres, each with a 5 km radius, placed 60 km off-center.To explore the parameter space, Lach et al. (2022) vary the central density, metallicity, and ignition offset, as well as add a few rigidly rotating models and one with a C-depleted core.Lach et al. (2022) find the familiar final state of a remnant surrounded by ejecta.The majority of the system, 1.09 -1.38 M ⊙ , remains bound.Lach et al. (2022) suggest their findings could account for Type Iax events of interme-diate and high brightness: they report ejected masses between 0.014 and 0.301 M ⊙ and ejected 56 Ni masses between 0.0058 -0.092 M ⊙ .They find that due to the density gradient inside the star, the deflagration only propagates outward from the ignition point, creating an asymmetry that gives a kickback velocity to the remnant.Depending on the magnitude of this velocity, the remnant may be ejected from the binary system.
Leung & Nomoto (2020) perform 2D studies of a C/O WD (and also C/O/Ne WDs, described in the following section) undergoing a pure deflagration using their own code Leung et al. (2015).This code uses a 5 th -order weighted essentially non-oscillatory (WENO) hydrodynamics solver, a level-set method to track the flame front coupled to a subgrid turbulence model, a three-stage burning model, a multipole gravity solver, and a Helmholtz EOS with Coulomb corrections.A 19-isotope reaction network is used to calculate the nuclear energy release.They simulate a quarter of the star, which has a cylindrical plus reflection symmetry.
Leung & Nomoto (2020) use a traditional 1.37 M ⊙ C/O WD with a central density of 3 × 10 9 g/cm 3 , a central temperature of 9 × 10 9 K, and solar metallicity which is added to the model as a 22 Ne mass fraction of 0.025.It is centrally ignited with an initially burned region with three perturbations, which are called "fingers".They explore the effect of changing the central density of the WD (0.50 -9.00 × 10 9 g/cm 3 ), the C/O ratio (0.3, 0.6, and 1.0), and the initial ignition type (central 3-finger and off-center single bubble).Their simulations eject a large amount of material (though the realization with the smallest central density ejected no mass); 0.92 -1.36 M ⊙ is ejected , 0.23 -0.34 M ⊙ of which is 56 Ni.The ejecta velocity is on the order of 5000 -8000 km/s.These simulations could account for the brightest Type Iax members, including 2005hk and 2012Z.

C/O/Ne WDs
There are not many studies that use these newly hypothesized hybrid progenitors, and the only 3D study available to date is that of Kromer et al. (2015), which uses an unmixed hybrid.Leung & Nomoto (2020)'s 2D study uses a mixed hybrid.Kromer et al. (2015) performed 3D simulations of deflagrations of a parameterized unmixed hybrid progenitor using an earlier version of LEAFS.They use a monopole gravity solver, and do not explicitly state the equation of state.
Their progenitor is a parameterized unmixed hybrid WD with a 0.2 M ⊙ 50/50 C/O core surrounded by a 1.1 M ⊙ , 50/47/3 16 O/ 20 Ne/ 12 C shell.A 0.1 M ⊙ accreted layer of 50/50 C/O is placed on top of this.The central density is 2.9 ×10 9 g/cm 3 .The resulting star is ignited off-center with an asymmetric 5 kernel ignition, and burning was suppressed in the O/Ne region.We emphasize that this parameterized model neglects the mixing that occurs during both cooling and simmering phases, which homogenizes the composition of the WD, typically up to 1.0 M ⊙ , before flame ignition.Also the suppression of burning in the outer shell ought to be expected to produce results different from the natural mixed state.
Their simulation results in a 1.39 M ⊙ bound remnant with 0.014 M ⊙ of ejecta.The ejecta is comprised of mostly unburned material, with 6.4 ×10 −3 M ⊙ of IGE, roughly half of which is 56 Ni.It has a low ejecta energy of 1.8 ×10 48 ergs of kinetic energy.This simulation gives a smaller 56 Ni-to-mass ratio than expected by observations, and also too fast rise and decline of the light curves.However, its ejected 56 Ni yields fall within range of that of the dimmest Type Iax events.
Unlike that of Kromer et al. (2015), the hybrid C/O/Ne WD used in Leung & Nomoto (2020) is mixed.It has a homogeneous composition of 20% 12 C, 50% 16 O, and 30% 20 Ne by mass, and then adds in solar metallicity as 2.5% 22 Ne.It has the same central density (3 × 10 9 g/cm 3 ) as their C/O model.They investigate the effect of changing the central density of the WD (1.00 -9.00 × 10 9 g/cm 3 ) and the initial ignition type (central 3-finger and off-center single bubble).They again find ejected yields similar to that of the brightest Type Iax events: ejecta masses in the range of 0.96 -1.30M ⊙ , 0.18 -0.36 M ⊙ of which is 56 Ni.These values are on average slightly smaller than their C/O WD counterparts.

This Study
This study seeks to determine if hybrid C/O/Ne WDs that mix during cooling undergoing deflagration in the single-degenerate, centrally ignited, near-Chandrasekhar mass scenario can reproduce the estimated yields of the dimmest Type Iax events.The main differences from previous works are the use of stellar evolution calculations to create the progenitor (Section 2.1) rather than a simpler parameterized model, and use of an ADR model flame in FLASH (Section 2.2), which is significantly different from the level-set method implemented in LEAFS and the software of Leung et al. (2015).We note that an ADR scheme in FLASH has been applied by several groups of researchers.The version we implement has had considerable development (see Section 2.2.2 and Appendix A.1) since an earlier implementation was applied to deflagration models of C/O WDs by Long et al. (2014).The results presented here are also directly compared to the DDT model from Augustine et al. (2019), which uses our same initial hybrid progenitor and software, as well as to observational data.The effect of the initial ignition region on the final yields is also explored.

METHODOLOGY
There are two main pieces needed for this study: a model of a hybrid C/O/Ne WD pro-genitor and software to perform the deflagration.Section 2.1 describes how the progenitor is formed from stellar evolution calculations using the software instrument Modules for Experiments in Stellar Astrophysics (MESA), and subsequently mapped to a 2D grid for use in simulation.Section 2.2 describes the software system, FLASH, used to simulate the deflagration.This includes a description of the most vital pieces: the equation of state (EOS) (Section 2.2.1); a module to model nuclear burning and the spread of the flame (Section 2.2.2); the initial conditions for ignition (Section 2.2.3); the initial parameters for the "fluff" material used in lieu of vacuum (Section 2.2.4); and the criteria for refinement of the grid (Section 2.2.5).In Section 2.3, we describe how a sample of 30 simulations (each called a "realization") was chosen to investigate the range of final yields, and how a subset of realizations was selected to explore the effect of the initial flame size on these yields.

Creating a Hybrid Progenitor
We used the same hybrid C/O/Ne WD profile as Augustine et al. (2019).This progenitor was created using the 1D stellar evolution software instrument Modules for Experiments in Stellar Astrophysics (MESA), version 8118, (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015)), by the following method, which evolves a zero age main sequence (ZAMS) star through the WD phase, followed by a period of cooling and accretion.A ZAMS star with a metallicity of Z=0.02 was evolved using the set of controls from Farmer et al. (2015) and the 29 species network "sagb NeNa MgAl.net".This network notably includes the 12 C( 12 C,p) 23 Na reaction of the carbon-burning process, which has been shown by De Gerónimo et al. (2022) to have a strong impact on the 20 Ne mass fraction of the end state WD.The ZAMS star evolves along the main sequence until it uses up the hydrogen in its core, transitions to the red giant branch (RGB), burns through its core helium, becomes a super-asymptotic giant branch (SAGB) star, and then ignites carbon burning off-center from the C/O core.This off-center ignition produces a carbon flame that is quenched by convective boundary mixing before it reaches the star's center (Denissenkov et al. 2013), leaving a core of unburned C/O fuel surrounded by a mantle of burned O/Na/Ne.The resulting configuration is a 1.09 M ⊙ hybrid C/O/Ne WD, with a 0.4 M ⊙ C/O core with a C fraction of 0.345, enclosed by a 0.69 M ⊙ O/Ne/Na mantle.
There is a Y e gradient between mantle and core (burned and unburned material), which in this model is ∆Y e = -0.00235(Brooks et al. 2017).The burned material has a smaller Y e due to the beta decay of 13 N that occurs during the burning, which motivates the choice of the extended "sagb NeNa MgAl.net" network.As described in Section 1.3, it has been demonstrated that this gradient causes the mantle and core to mix as the WD cools, reducing the central C fraction and homogenizing the WD's interior far before accretion and subsequent ignition would occur (Brooks et al. 2017;Schwab & Garaud 2019).
This mixing significantly decreases the carbon fraction in the core from 0.345 to 0.1419.
To model the accretion phase, the WD was placed in a binary system with a 1.4 M ⊙ He star with a 3 hour orbital period, following the prescription of Brooks et al. (2016).The He star was also evolved from a ZAMS star with a metallicity of Z=0.02, except it used the 8 species network "basic.net".The accreted layer therefore is devoid of the heavier elements, but this isn't expected to change the deflagration dynamics, which occur in the center of the star at the temperature peak.
As the WD gains mass through accretion from its He-rich companion, it compresses and increases in temperature.The mixed central part of the WD becomes hot enough to start "simmering": C-burning activates, but convection is able to carry excess heat throughout the material and prevent thermonuclear runaway.The simmering decreases the C-fraction a bit further, to ∼ 0.139.Beyond that, the underlying accreted layers begin burning He, enriching the outer layers of the star with C and O. Finally, the newly accreted material rests on the surface.Accretion continued until the convection in the core could not disperse the excess heat fast enough and thermonuclear runaway was imminent; in MESA this was considered to be the point where the central temperature was on the order of 10 8.9 K.The resulting star just before the deflagration has a mass of 1.36 M ⊙ , with a central density of 2.2 × 10 9 g/cm 3 and a central temperature of 8.2 × 10 8 K.
The resulting progenitor WD was mapped to a cylindrically symmetric 2D grid to use as the input for a deflagration simulation in FLASH.The code that performed this transformation can be found on github (see Software section), and the transformation is described in Willcox et al. (2016) and was performed for Augustine et al. (2019).
First, the profile was mapped to a 1D, 4 km grid (the maximum resolution of the simulation, as described in Section 2.2).Next, the 29 species in the MESA profile were condensed into 4 in order to be used with our burning model in FLASH: 12 C, 16 O, 20 Ne, and 22 Ne.The prescription for doing this is to keep the 12 C mass fraction (X12 C ) the same, and then use 22 Ne as a proxy for the neutron-rich nuclides by calculating its mass fraction (X22 Ne ) using Y e .The remaining mass fractions of 16 O (X16 O ) and 20 Ne (X20 Ne ) were kept in the same ratio as in the MESA profile.
To force the grid to maintain hydrostatic equilibrium, the cell densities were adjusted.The density, temperature, and composition profile of the mapped model can be seen in Figure 1. and composition (bottom) of the end state of the hybrid C/O/Ne WD after accreting mass from its companion star.The 29 element reaction network was condensed into 4, so it can be used with our burning model in FLASH.Here, 22 Ne is used as a proxy for neutron-rich material, and the accreted He has been added as 20 Ne as they have the same electron-to-baryon ratio.
Three "layers" can be seen, with transitions between them.The first is the compositionhomogenized core convection zone (1.07 M ⊙ , 1.018×10 8 cm), which has enveloped the cooled, mixed star and further lowered its C-fraction.The expected adiabatic temperature profile can be seen in this layer.The next layer extends out to 2.15 ×10 8 cm, and contains ashes from He-burning, where the accreted material has compressed and become hot enough (1.5-2.5 ×10 8 K) to burn He.The remaining 10 −6 M ⊙ layer that extends to the edge of the star (2.462 ×10 8 cm) is made of unburned accreted material.Perhaps the most obvious difference between the MESA and FLASH models is the large amount of 20 Ne in the accreted layerin reality, the accreted layer is mostly 4 He, but this is not an element in the burning model.Since 20 Ne and 4 He have the same electron-tobaryon ratio, this should not change the equilibrium much.

Simulating a Deflagration
The deflagration simulation was performed using a customized version (Calder et al. 2007;Townsley et al. 2007Townsley et al. , 2009Townsley et al. , 2016;;Willcox et al. 2016) of FLASH, a modular, extensible, highly parallelizable software system (Fryxell et al. 2000;Calder et al. 2002).FLASH allows the combination of physics modules to create unique problems, which are then solved on a grid.
For this problem, we use FLASH version 4.6.2 with the PARAMESH adaptive mesh grid with a maximum resolution of 4 km; the split PPM hydro solver (Colella & Glaz 1985); the Poisson multipole gravity solver; a Helmholtz equation of state (EOS) with Coulomb corrections that describes the degenerate electron-positron plasma of the WD, overviewed in Section 2.2.1; and a burning model specifically designed for modeling the thermonuclear processes in Type Ia supernovae that is described in Section 2.2.2.For compilation, we used the GCC 10.2.0 compiler, MVAPICH version 2.2rc1 compiled with GCC, and parallel hdf5 version 1.8.20.

Equation of State
The EOS used for our application describes the degenerate electron-positron plasma that makes up the WD.It is the Helmholtz EOS described in Timmes & Swesty (2000) with addi-tional Coulomb corrections.Details of the implementation of this EOS are collected in Appendix A.2 and described in brief in Fryxell et al. (2000).
As a deflagration produces slow-moving, sparse, and cool ejecta, it is important to be aware of the limits of this EOS.The Helmholtz EOS is meant to describe a fully ionized plasma, and its range of use has been cited as 10 −10 < ρ < 10 11 g/cm 3 and 10 4 < T < 10 11 K (Timmes & Swesty 2000).However, there is a region in ρ -T space in which the Helmholtz EOS will return a negative internal energy.The bounds of this region are the right-triangle-shaped contour E int ∼ 10 11 ergs, which lies where T <∼ 10 6 K and ρ <∼ 10 3 g/cm 3 .This region extends to lower densities as the temperature drops.This problematic region is also present in the newest improvement of the Helmholtz EOS, the Skye EOS, which also describes a fully-ionized plasma over an even wider temperature and density range (Jermyn et al. 2021).Jermyn et al. (2021) find a triangular region of parameter space that fails when calculating the thermodynamic consistency relations, which they concisely explain: "The feature at intermediate densities and low temperatures indicates negative pressures caused by our assumption of a fully ionized free energy in a region that should form bound states, indicating that Skye is not valid in that limit."We discuss workarounds to this issue in Section 2.2.4.

Burning Model
The burning model has been developed from its original version in Calder et al. (2007); Townsley et al. (2007) to the current implementation of Townsley et al. (2016), and continues today.It is based upon the advectiondiffusion-reaction (ADR) scheme of Khokhlov (1995) and Vladimirova et al. (2006).The addition of 20 Ne to the model, which is necessary to simulate burning of this hybrid progenitor, is described in Willcox et al. (2016).Summarized below are the key characteristics, and more detail can be found in Appendix A.1 and the above references.
A 3-stage burning process is used to model the evolution and energy release of a physical flame: carbon burning, burning to Si-group dominated nuclear statistical quasi-equilibrium (NSQE), and burning to Fe-group dominated nuclear statistical equilibrium (NSE).Each of these three burning stages is tracked by its own scalar progress variable -ϕ f a (fuel to ash), ϕ aq (ash to NSQE), and ϕ qn (NSQE to NSE), respectively -that represents the fraction of a grid cell that has undergone that type of burning.It logically follows that 1 ≥ ϕ f a ≥ ϕ aq ≥ ϕ qn .ϕ f a is evolved using an ADR scheme, and ϕ aq and ϕ qn at tuned timescales (Calder et al. 2007;Townsley et al. 2016).The properties of the final NSE state are evolved to account for changes in nuclear binding energy due to electron-capture as well as changes in composition due to hydrodynamic evolution.Willcox et al. (2016) found that the addition of 20 Ne did not require any substantial changes to to the structure of this model. 20Ne is consumed in the first (C-burning) stage, and only serves to make the burning slower and less complete.The remainder of the stages progress with the same timescales as previously determined.Therefore, the only required change was to add 20 Ne as an element in the composition of fuel for the first stage of the burning model.
Since the length scales of a realistic flame cannot be simulated inside a WD, the model flame is artificially thickened so that the flame front is exactly 4 grid zones wide.The minimum resolution needed for an accurate model is 4 km (Townsley et al. 2009).The laminar flame speed is interpolated from a table of density, 12 C-fraction, and 22 Ne (a proxy for metallicity) fraction, and then boosted to account for unresolved Rayleigh-Taylor instability.To obtain laminar flame speeds for the lower densities and C-fractions of our hybrid progenitor, the table is simply extrapolated (Willcox et al. 2016).

Ignition
As the actual flame cannot be resolved in this simulation, the ignition is approximated by a fully burned region whose surface is perturbed to capture the small-scale effects of convection.The deflagration was initialized at the temperature peak of the progenitor, at the center, with a 150 km radius "matchhead" of fully burned material (ϕ qn = 1).The shell of the matchhead was perturbed by an angular sinusoid function with a random harmonic number and amplitude, following the prescription of Townsley et al. (2009): The values of ℓ are taken from Townsley et al. (2009): ℓ max was chosen so that that the perturbations are resolved, and a relatively large choice for ℓ min allows the star time to expand and more accurately reproduce DDT observational yields.The solutions are limited to those that have only inward radial perturbations at the cylindrical symmetry axis, to prevent long thin plumes from forming and growing along the axis, a nonphysical behavior caused by our choice of cylindrical coordinates.In 2D simulations such as this, m = 0 only is dictated by the assumed cylindrical symmetry.The sum reduces to: These A ℓ are set at the beginning of the simulation by a random seed.As will be described in Section 2.3, for a statistical approach we run a suite of simulations, each of which has a different random seed that generates a unique set of A ℓ .

Initial parameters
As the simulation is in 2D, a simplification that readily follows is to start with a zero background velocity field.In reality, the core of the star is convective, and future 3D simulations will include a velocity field derived from turbulence simulations (Jackson et al. 2014), but as 2D simulations cannot describe turbulence, we neglect an initial velocity.
As the PPM hydrodynamics solver of FLASH requires material on every cell of the grid, the vacuum of space is instead replaced with a lowdensity "fluff" material surrounding the star (Zingale et al. 2002).The fluff composition was set to match the abundances of the outermost zone of the initial profile, to eliminate any composition gradient.The fluff pressure was set to 5×10 14 dyne/cm 2 , 25% higher than the pressure at the edge of the star.
The pressure of the outermost zone of the initial stellar profile is resolution-dependent.The steep edge of the surface of the star results in a resolution dependence of the surface pressure.We found that within the first timesteps of the simulation that the pressure at the edge of the star can jump up an order of magnitude from ∼4 × 10 14 dyne/cm 2 to ∼4 × 10 15 dyne/cm 2 , which necessitated setting the fluff pressure to mitigate any resulting outwards pressure waves and subsequent inward shock due to this change.Since our deflagration is slow moving and the ejecta densities are sparse, we used a fluff density of 10 −4 g/cm 3 .
During the expansion of the ejecta after the explosion, some regions pass through the problematic region of our equation of state discussed in Section 2.2.1 near a density of 10 −2 g/cm 3 and temperature of 10 6 K. Interestingly, the areas where this negative internal energy are returned are the grid cells immediately in front of the plumes created by the flame, which have a much lower temperature (∼10 6 K) than the surrounding material.It is important to note that while this may be physical, it is likely due to an artifact of the piecewise parabolic split hydrodynamics solver of FLASH.The cause of these low temperature regions will have to be further investigated.
In the meantime, to attempt to avoid this area of parameter space, we set the relatively high temperature floor of 2 × 10 6 K.This high temperature floor introduces a negligible amount of internal energy into our simulations at later times (∼ 10 s).Since burning has stopped at this point and internal energy is not included in the bound vs unbound criteria, this increase in temperature floor does not significantly change the results.This is clearly an imperfect solution and some simulations were not able to successfully continue to the desired stopping time and were ended early.

Refinement
We have chosen a fluff density on the order of 10 −4 g/cm 3 and the ejecta reaches densities of the order 10 −2 g/cm 3 or lower.As a result, distinguishing between fluff and star material based on density, as was done in Townsley et al. (2009), becomes a less tenable strategy toward the end of our simulation.Instead, we replaced this with a combination of a fluff marker and a derefinement scheme.Adding a fluff marker was practically accomplished by initializing the simulation with a scalar variable, 0 ≤ ϕ fluff ≤ 1, where ϕ fluff = 0 inside the star, and ϕ fluff = 1 outside.This scalar is advected by the hydrodynamics much like a single species abundance, so that it effectively tracks the percentage of fluff in each grid cell.Rather than principally using density to detect non-star material, this marker was used to keep the fluff as derefined as possible, (neighboring blocks can only differ by 1 level of refinement at most), by setting the maximum refinement level of the fluff to 1 (4096 km), the largest grid cell size available for our simulation.
Burning stops at about 3 s of simulation time due to the decreasing density leading to the extinction of the propagating flame.The flame and dynamic ash were manually turned off at 5 s (some simulations burn faster or slower than others, so 5 s was chosen so as to not prematurely cut off the burning).The flame needed to be turned off because, as discussed in Section 4.9, the bound remnant in some simulations contracted to high enough densities that burning might occur.We have opted to consider such cases separately since the necessary conditions for reignition and the nature of the combustion under those conditions are uncertain.The densities and temperatures reached by the material in our simulation are low enough where continuing with dynamic ash enabled makes a negligible difference.
At the start of the simulation, the maximum refinement level was set to 11, which corresponds to a grid cell size of 4 km, the size necessary to accurately model the flame.The maximum refinement level of the non-energy generating regions (ie anything not burning) was set to 9 (16 km), and for fluff 1 (4096 km).After 3 s, for derefinement, we decreased the non-energy generating regions' maximum refinement level by 1 (ie doubled the minimum grid cell size) every time the star doubled in size: 3 s (refinement level 8, 32 km), 5 s (refinement level 7, 64 km), and 10 s (refinement level 6, 128 km).Since after 3s burning has stopped, all regions are classified as non-energy-generating and are therefore subject to this refinement constraint.

Suite
To obtain results representative of a population, rather than a single member, we ran a suite of simulations with different initial flame perturbations, each of which we call a "realization" (see Section 2.2.3).It has been demonstrated by Calder et al. (2010) that 15-20 realizations of a DDT simulation are enough to cause the standard deviation of the 56 Ni yield of the sample to converge, however it is unclear whether this applies to deflagration simulations.Nevertheless, to follow the prescription of Augustine et al. (2019) which used the same hybrid progenitor undergoing DDT, thirty realizations (numbered R11-R40) were performed for this study.We believe that this provides useful, though incomplete, sampling of possible outcomes.
To create a more direct comparison to WD undergoing DDT, each realization used the same random seed as its DDT counterpart in Augustine et al. ( 2019), so that the geometry of the initial flames were perfectly replicated.The initial conditions of the fluff and the refinement criteria, however, had to be adjusted from that of Augustine et al. (2019) as described above in Section 2.2.4 so that the deflagration could evolve for the requisite longer period of time.A description of how the data from these thirty realizations was analyzed can be found below in Section 3, and a summary and comparison of the results of these thirty realizations is shown in the following Section 4.
To better understand the effect of the initial flame radius on the results, we performed additional simulations using 10 realizations from our sample that represented the full range of yields.The flame radius was changed from its original 150 km to 75 km ("half-radius") and 300 km ("double-radius").

ANALYSIS
Analysis was performed and plots were created using the yt package (Turk et al. 2010).For each realization, data was collected every 0.1 s of simulation time.To distinguish between the remnant and ejecta, we separate our results into regions of positive and negative gravitational binding energy, respectively, calculated by Equation 4. We choose this separation so that we can directly compare these simulated yields with those estimated from observational data, as discussed in Section 4.8.
Here, E grav is the gravitational potential energy, and E kinetic is the kinetic energy (KE) of each grid cell.Internal energy was not included in this calculation because the high temperature floor artificially increased the internal energy of the low density, low temperature regions at the edge of the plumes of burned material, which are often located near the line of 0 binding energy.The internal energy of these regions is orders of magnitude smaller than the kinetic or gravitational potential energy at late times, so the final yields of the remnant and ejecta are negligibly affected by excluding internal energy from the binding energy calculation.
Since our combustion model does not explicitly evolve abundances of particular nuclear species, as discussed in Section 2.2, the yields of specific elements, namely IGEs and 56 Ni, must be calculated using the progress variables (ϕ f a , ϕ aq , ϕ qn ), neutron excess of the fluid (Y e ), and initial, unburned mass fractions of 12 C (X12 C ), 20 Ne (X20 Ne ), and 22 Ne (X22 Ne ) the proxy for neutron-rich material).
The simulation considers the remaining material as 16 O.
The mass of material in each grid cell was calculated using the mass fraction, density, and 3D volume of each cell.
where M i is the mass of substance i, X i is the mass fraction of substance i, ρ is the density and V the 3D annular volume of a cell.The mass fraction of IGE is simply the progress variable ϕ qn , as described in Section 2.2.2, since the difference between baryon density and mass density is unsubstantial.
The calculation of the mass fraction of 56 Ni present in each grid cell (X56 Ni ) follows the estimation of Equation 29 of Townsley et al. (2009), which assumes that 54 Fe and 58 Ni are the next most-common elements in the fully burned material.where Y e,i is the electron fraction of substance i.The electron fraction in the fully burned material is and the electron fraction in the unburned fuel is

SIMULATION AND RESULTS
The evolution of the simulation is described in general in Section 4.1, with explanations of the most noteworthy features and their connection to the physical system.This is followed by a description of how the final system is split into two parts -remnant and ejecta (Section 4.2).The final yields are tabulated and discussed (Section 4.3) alongside notable correlations (Section 4.5) and an analysis of the changing kinetic energy of the system (Section 4.6).These results are then compared to a DDT simulation by Augustine et al. ( 2019) that uses the same progenitor (Section 4.7), as well as to observational data (Section 4.8).

Overview of Explosion Stages
Figure 2 shows both the density and percentage of material burned to IGE (ϕ qn ) of an average realization, R25, at 0, 0.7, 3, 13, and 20s of simulation time.Within the first second, the surface of the star puffs up and mixes with the fluff, creating artistic but entirely unphysical tendrils of swirling material.This effect is illustrated in the last three panels of Figure 2. The initial matchhead is put into pressure equilibrium with the surrounding star, but it is still a perturbation to the initial HSE as it introduces an early change in composition and density.The outward propagation of sound waves due to this perturbation is what we believe causes this sudden surface expansion.
In the meantime at the star's center, the ADR scheme described in Section 2.2.2 starts to propagate the model flame, approximating the physical flame which heats the surrounding cold fuel by conduction until it begins to fuse.The resulting released nuclear energy is absorbed by the newly burned material as heat, increasing its temperature and causing it to expand, which decreases its density.This leaves pockets of hot, less dense material surrounded by cold, denser fuel.The Rayleigh-Taylor (RT) instability at this boundary perturbs the flame surface, causing plumes to form, merge, and wrinkle, which increases the surface area of the growing flame so it spreads faster.Our model captures the geometry of this RT instability at a larger resolved scale (∆x = 4 km), and adjusts the flame speed accordingly to account for unresolved small-scale motion.This strategy has known limitations due to the absence of realistic turbulence in 2D.
The plumes begin to rise and accelerate due to buoyancy, pushing the unburned fuel out of the way and imparting it with kinetic energy.The burned region is also increasing in size due to the flame propagating.The effect of these processes is an overall expansion of the star.As the flame front propagates and plumes rise, the density of the unburned fuel the flame propagates into decreases.A physical flame in a lower density environment would become slower and wider, since it takes longer to burn material at lower densities.At some point, the corresponding steady-state flame at a density of around 5 × 10 7 g/cm 3 would grow wider than the radius of the star.Our artificial flame is always 4 grid zones wide, so these effects are modeled by a decrease in the laminar flame speed and an increase in the evolution timescales of the progress variables ϕ aq and ϕ qn .Eventually the timescales become so long that material essentially stops burning, a bit before 3s.After this point, there is no more energy being added to the system from fusion, and the expanding material becomes ballistic.Gravity will slow this expansion.Unlike DDT models, where the supersonic detonation produces a fairly uniform and rapid expansion, the system remains dynamic, and there are motions and interactions between the burned and unburned material.
The low carbon fraction in the star's center means that less binding energy is released during fusion than in a regular C/O WD.This causes the flame to propagate more slowly, allowing the star more time to expand and leading to less buoyant plumes and more incomplete burning.Most of the star remains unburned, only now more mixed and less dense due to expansion.Because of this, most of the star remains behind and its gravitational pull causes most of the surrounding material to fall back, creating a remnant of unburned fuel mixed with ashes from burning.This fallback causes small shockwaves that push mass back out, but do not give it enough energy to fully escape the gravitational force, so that the mass falls back again and the process repeats.This behavior will hereafter be referred to as the pulsation of the remnant.Simultaneously, in contrast to this fallback, the puffed up surface of the star as well as the outer layers of the burned plumes have enough kinetic energy to escape the gravitational pull of the remnant, and travel out into space, continuing to decrease in both temperature and density.

Stopping with Separated Ejecta and Remnant
The system is separated into two pieces -remnant and ejecta.It is assumed that all the bound material will eventually fall back and become part of the remnant, while the ejecta will move off into space to be detected by our instruments on Earth.Therefore, the ejecta yields alone will be compared to yields calculated from observational data.We follow the evolution of this remnant-ejecta system out to 20s of simulation time, when the ejecta mass has become constant with time.Too much farther in time than this point, the ejecta and the fluff would begin to interact, and the fluff would decelerate the ejecta, imparting a nonphysical effect that would give an underestimate for the ejecta kinetic energy.This could be mitigated by implementing a density rolloff of the fluff.Without a detonation to shock the surface layers of the star, this becomes the biggest challenge to continuing this simulation to longer times.
Another issue with running these simulations out to long times is the disparate time evolution between the ejecta and the remnant.The ejecta is relatively cool, sparse, and settling into free expansion in a fixed form.In contrast, the remnant is hot and dense, collapsing in on itself and mixing on short timescales.This behavior is apparent in Figure 2; the expanding plumes remain essentially constant in shape from 13 s to 20 s, while the remnant is gaining mass and starting to pulsate.Resolving the simulation at the timescale of the remnant is a severe limitation on the evolution of the ejecta.As our primary interest for this 2D study is in the ejecta yields, 20 s is long enough to obtain the desired information, but this issue will require additional care when running future 3D simu-lations and creating lightcurves and spectra, as discussed in Section 6.
Other studies (Kromer et al. 2015;Lach et al. 2022) deresolve the remnant in order to observe the ejecta at later times, but even if we followed this prescription, our simulations would still be forced to stop short.While the majority of our realizations reach the full 20 s, some run into the region of parameter space described in Section 2.2.4 that our EOS is not made to handle, even with the increased temperature floor.As we did not want to further increase the temperature floor and put more energy into our simulation, and as the ejected masses had leveled off already, we did not run these realizations further.All simulations made it past 15 s.The data reported in Section 4.3 are different quantities at the end of our simulation.It is important to note that these final yields refer to the yields at the end of each realization, which is either 20 s or at the point where the simulation stopped due to the above described EOS difficulty.The accompanying correlations between variables should not be affected by the difference in runtimes however, as the quantities of interest have leveled off, and should not change much over time.

Final Yields
The ranges of ejecta and remnant yields and energies are shown in Table 1, visually supplemented by the box-and-whisker in Figure 3.The table reports the full range and also the middle 50% (the interquartile range, or IQR) of the data.Each box-and-whisker plot shows the median (50th percentile, orange line) of each quantity, with upper and lower box boundaries marking their 75th and 25th percentiles, respectively.The range enclosed by the box corresponds to the IQR of the data, which is also reported in Table 1.Note that the median of the data is not necessarily equidistant from the 25th and 75th percentiles (ie not necessarily in the center of the box).The whiskers are plotted  4) having a negative or positive value, respectively.The IGE, NSQE, Ash, and Fuel yields correspond to X N , X q , X a , and X f described in Appendix A.1 and Section 2.2.2, and the estimate of 56 Ni was calculated using Equation 8.It is important to note that the minimum values listed do not all come from the same realization (and the same with the maximum values).
at 1.5 × IQR away from the 75th and 25th percentiles, but do not extend beyond the limits of the data.Data that lie outside the range of the whiskers are plotted as individual points.Notably, ∼ 66 -77% of the star remains unburned; ∼ 14 -21% is fully burned to IGE; only ∼ 2 -4% is partially burned to NSQE; and ∼ 7 -12% undergoes carbon burning and then stops (the ash state).Only ∼ 5 -9% of the star is burned to 56 Ni, which is ∼ 40% of the IGE.In addition, the vast majority, ∼ 90 -99 %, of the star remains bound, demonstrating that a deflagration does indeed produce much smaller amounts of ejected material and high-mass el-ements than a DDT.The small amount of remaining mass has enough energy to overcome the gravitational pull of the remnant, and is ejected from the system.As demonstrated in Figure 3, most of this ejected material is unburned fuel, but it also contains the edges of the plumes of burned material.Interestingly, the majority of the burned material, including the IGE and 56 Ni, is bound, and will fall back and join the remnant.
One important feature to note is that the simulation that produced the largest ejected mass (R16) is neither the one corresponding to the most ejected 56 Ni mass (R20), nor the one with the most ejected IGE (R14), nor the one that burned the most mass (R19).Similarly, the simulation that produced the least ejected 56 Ni and ejected IGE (R32) has neither the least ejected mass (R35) nor the least burned mass (R28).This feature is more clearly demonstrated in Figure 4, which shows how these quantities evolve with time.
In these plots, the results of the following six simulations are highlighted: R16, which produced the most ejected mass; R20 and R32, which yielded the most and least ejected 56 Ni mass, respectively; R28, which by far burned the least mass; and R23 and R25, two average runs that produced almost the same ejected 56 Ni mass but exhibit different rise times.Highlighted in each plot is the full range of the data from the suite, and also shown is the region where 75% of the data is distributed, and the region where 50% of the data is distributed (the IQR).These measures were chosen rather than the standard deviation, because they give a better picture of how our realizations are distributed.The results are not perfectly Gaussian, and have a large range; therefore the standard deviation is often greater than the average, and shading in that region would give a misleading representation of our data.Therefore, we use quartiles as a better description of our data distribution.It can be seen that the full range of the data is about twice as large as the middle 50% of the data.
In general, the more mass that is burned by the flame, the more material that is converted to IGE, and the more nuclear energy that is released, which leads to higher kinetic energies imparted to the surrounding star, and therefore more material that can break free from the remnant and is ejected.While this reasoning seems straightforward, things are complicated by both the nonlinearity of this problem, especially the RT instability, and the overly restrictive cylindrically-symmetric 2D geometry.The RT instability determines how the initial perturbations evolve into plumes of burned material, which plumes overtake others and which become dominant, and how the plumes combine, deciding what parts of the star will contain the most burned material.It also decides the initial direction that these plumes move, which means they are either headed towards the less dense, colder material towards the edge of the star, or towards the hot dense fuel of the inner star, which decides how long the plumes can burn.
Our choice of geometry, on the other hand, causes the plumes to move away from the z-axis and curl around the mostly unburned core of the star.In our cylindrically symmetric geometry, each grid cell corresponds to an annulus of material, which means that grid cells along the r-axis occupy a larger volume and contain more mass than those along the z-axis.This means that plumes parallel to the z-axis tend to rise faster and burn through less mass, since they have less material overhead.This creates a shear vortex that rolls the plume heads towards the r-axis (Townsley et al. 2009).As discussed in Section 2.2.3 the initial flame morphologies are explicitly chosen to not have a radial maximum on the z-axis axis because that direction has effectively too little resistance to plume rise due to the imposed symmetry -allowing a maximum there leads to excessively unrealistic pathologies.This curving is apparent in R28, shown in Figure 5, which has the least total burned mass and the least total IGE mass, but ejects an amount of mass and IGE within the middle 50% of all the realizations.

Sensitivity to Initially Burned Region
It is important to quantify the effect that the shape and size of the initially burned region (matchhead) has on the results, as the matchhead is only necessary because the realistic initial flame scale cannot be resolved.As seen in Section 2.2.3, the matchhead is described by Equation 3, where each realization has a different set of A ℓ .
Due to the nonlinearity of this problem, it is nearly impossible to relate the shape of the matchhead to the ejecta masses.The resulting ejecta mass is strongly correlated with neither the initial burned mass (a measure of the volume of mass inside the flame at t=0), nor the maximum or minimum radius the flame plumes reach.Figure 6 shows the realization with the least ejected 56 Ni (R32) side-by-side with the realization that produced the most (R20).The burned regions look similar at 3s, but R32's main plume is curling inward, whereas R20's plumes are all moving outward.
This difference shows that the early flame evolution can create large differences in outcome, but it is difficult to predict exactly which initial parameters and in what combination produce the effect.This lack of correlation between the initial flame geometry and the final yields shows that the choice of initial conditions for the simulation suite (discussed in Section 2.2.3) successfully creates a random sample in the sense that the shape of the initial flame cannot be directly related to the yields.
However, changing the radius of the initial matchhead does produce a systematic effect on the results.From our main suite of 30 realizations, we selected 10 that spanned the full range of ejected 56 Ni yields.We investigated how doubling (300 km) and halving (75 km) the matchhead radius impacted the final yields, using the same prescription as before in Section 3 to analyze the results.
Examination of this data yields multiple interesting results.A spatial comparison of the results of these two matchhead radii for the same realization (R20) can be seen in Figure 7. Perhaps the most prominent feature is that the larger matchead produces much smaller and more constrained plumes of burned material.This is because the larger matchhead deposits much more energy into the star's center at the start of the simulation, causing the star to expand faster.Additionally, the larger matchhead suppresses the flame instabilities that boost the burning; the smaller matchhead's plumes will evolve in the dense center of the star subject to the RT instability, while the larger matchhead has a flat profile in the denser region.Due to this expansion and suppression, a larger radius matchhead burns through less dense mate-rial more slowly, causing the burning to extinguish sooner and less material to be burned to IGE than for a smaller matchhead.This means the buoyant plumes of burned material will not only be smaller in size, but also will rise less because of the lower density of material surrounding them, leading to a smaller overall KE.This can be seen in Figure 7; the smaller matchhead's plumes reach to the edge of the expanding star,  4, the results of the following four simulations are highlighted: R20 (yellow) and R32 (black), which yielded the most and least ejected 56 Ni mass for the baseline (150 km) radius matchhead, respectively; and R23 (blue) and R25 (fuchsia), two average runs that produced almost the same ejected 56 Ni mass but exhibit different rise times for the baseline (150 km) radius matchhead.Also shown are the IQR for each timestep (purple, shaded), the range of the middle 75% of the data (blue, shaded), and the full range of the data (light blue, shaded).
while the larger matchhead's plumes are much more constrained.
As shown in Figure 8, as the initial matchhead radius increases, the total kinetic energy, 56 Ni, and IGE yields decrease, as does the 56 Ni:IGE ratio.Additionally, the burning stops sooner.The final yields from these two different matchhead radii are shown in Tables 2 and 3. Most  50% 2.624 -4.520 1.114 -1.515 1.552 -3.193Table 3.The range and interquartile range (IQR, middle 50%) of the final yields from the doubled (300 km) matchhead radius simulation suite.
is that eight of the ten double matchhead radius realizations produced no ejected mass.In these instances, all of the mass remained gravitationally bound and fell back towards the progenitor, suggesting that the white dwarf will re-form without any net mass loss.
It is important to note that all of these initial flames, regardless of size, are but an approximation, as a realistic flame is smaller than the size of a grid cell and therefore cannot be resolved.Changing the initial size of the matchhead can be thought of as starting the simulation at a different point in time; a realistic flame would need less time to grow to the size of the smaller matchhead than the larger matchhead.The larger matchhead also expands the star faster and has suppressed flame instabilities.Logically, a smaller initial flame seems like it would be more realistic, but it is uncertain whether the ignition would be in one location and grow, or from several locations and then merge.Other studies test the differences between this by performing a similar initial conditions study with multiple ignition points and sizes, but this is beyond the scope of this work.A flame in 2D like ours cannot be turbulent, so is not realistic.Therefore, this test cannot say anything about the "realism" of the initial parameters of the flame, but rather serves as a measure of uncertainty for our simulation.Though the initial matchhead size produces a systematic effect, the yields remain in the suggested range of Type Iax events.

Correlations
Returning now to our suite of 30 realizations, the nonlinearity and curvature of the simulation is one of many reasons why, even though there is the expected correlations between the burned mass, the ejected mass, and the IGE mass, they exhibit varying degrees of tightness and scatter.These correlations, among others, are shown in Figure 9.The burned mass is made up of ash, NSQE elements, and IGE which in-cludes 56 Ni; therefore it has a positive correlation which each of these quantities, as shown in the first column of Figure 9.These correlations display a scatter which generally increases with increases burned mass.Additionally, the more mass that is burned, the less total mass remains bound, though this is not a strong correlation.The IGE and 56 Ni show a strong positive correlation, meaning that a reasonably well-defined 56 Ni:IGE ratio can be calculated.This ratio is ∼ 0.4 for total, bound, and unbound quantities, but has a much larger range for the unbound quantities, as seen in Table 1.The perfect correlation between the total bound mass and the total unbound mass is due to the fact that the bound mass plus the unbound amss of the system must sum to the total mass, 1.36 M ⊙ .Unbound quantities tend to exhibit a stronger correlation than bound quantities, for example the graph of unbound 56 Ni and unbound IGE is more strongly correlated than that of bound 56 Ni and bound IGE.However, the unbound quantities tend to increase their scatter with increasing burned mass much more so than their bound counterparts.The quantity that exhibits the strongest determinant as to how much material will be bound vs unbound is the system's kinetic energy.

Kinetic Energy
The pulsation of the remnant can be seen in the maximum density vs time as well as the total KE vs time.The latter is plotted in Figure 10.
Both quantities exhibit a damped oscillatory pattern.The maximum density quickly drops two orders of magnitude during burning, as the energy deposited by fusion expands the burned regions and pushes the rest of the star out with it.The densest part of the simulation is still the expanded core of the star.In general, realizations with less burned mass will reach higher densities faster, because the star is less disrupted by the burning, less expanded, and can therefore recombine into a remnant faster.Realizations that burned less mass also impart less KE, so they reach a smaller peak KE sooner.This picture is complicated by the fact that the KE is not evenly distributed throughout the star, and some simulations are left with more bound mass than others.R35 reaches the lowest peak KE and it has the most bound mass, whereas R13 reaches the highest peak KE but R16 has the least bound mass.
As shown in Figure 10, the total KE reaches its maximum right around the time burning stops, since the energy released by fusion is what causes the temperature of the burned regions to increase, leading to expansion of the plumes that push mass out from the center of the star due to buoyancy.As soon as the burning stops, the total KE of the simulation starts to decrease, as it becomes converted into gravitational PE.
In general, more burned mass means more energy released and therefore a higher maximum KE.However, as shown in the bottom panel of Figure 11, the amount of burned mass alone is not enough to determine the maximum KE of each realization.While there is a correlation between these two quantities, there is scatter, indicating that there are other factors at play.Due to the nonlinearity of the simulation, realizations that have a certain amount of burned mass will impart different amounts of KE depending on the location and direction of the plumes.This is also why the final burned mass is not strongly correlated with the ejecta yields.After the burning has stopped, however, the nonlinear effects from the RT instability and curving due to the geometry are not as prominent.In general, the kinetic energy imparted during the burning phase determines how much mass will be able to escape the gravitational pull of the remnant.Therefore, the maximum total KE is tightly correlated with the final unbound mass, as seen in the top panel of Figure 11.There appears to be a "characteristic" ejecta velocity, since the ejecta KE and ejected mass show a linear relationship, as seen in Figure 12.Using the slope of the line of best fit, this characteristic ejecta velocity can be estimated as 3500 km/s.

Comparison to DDT model
The initial white dwarf progenitor model, initial flames, burning model, and simulation environment (FLASH) are identical to those in Augustine et al. (2019)'s DDT study, while the initial conditions, refinement criteria, and explosion mechanism differ.These differences might not allow for a detailed realization-byrealization comparison, however a general comparison between the two explosion mechanisms can certainly be made.
At the beginning of both simulations, the burning rate and IGE production closely follow each other, as the progenitor is undergoing a deflagration.The deflagration transitions into a detonation at about 1-2s of simulation time; realizations that burn less mass during the deflagration phase tend to transition sooner.
As a detonation incinerates the entire star, it is expected that our WD progenitor undergoing a DDT will have higher burned yields and explosion energies than when undergoing a deflagration alone.This large difference is shown in Figure 13, which shows the disjoint ranges of final total IGE and 56 Ni mass produced by the two explosion mechanisms.The DDT simulation, besides producing over 4 times greater yields, also has a much larger range of yields than its deflagration counterpart.Additionally, the 56 Ni:IGE ratio of the DDT simulation is around 0.75, nearly twice than that of our deflagration simulation, 0.4.

Comparison to Observations
Using Arnett's model (Arnett 1982) or Valenti et al. (2007)'s variation of it for a inner highdensity region surrounded by a lower-density region, ejected masses can be calculated using parameters inferred from observations, including the ejecta velocity, opacity, and maximum brightness of the event.Table 4 collects the results of these calculations performed by other studies for the dimmest Type Iax events.Where possible, data from multiple studies of the same Type Iax event is included.It is important to note that while there are over 70 members of the Type Iax subclass, there are only 13 events listed in Table 4.This is because generating lightcurves and spectra are enough information to classify an event as Type Iax, and analysis to determine the ejected masses and energies requires several additional assumptions and calculations.This search through the literature is not necessarily complete, but serves to demonstrates the range of brightnesses of Type Iax events: 0.001 -0.30M ⊙ ejected 56 Ni.
In our simulation, it is assumed that all the bound material will eventually fall back and become part of the remnant, while the ejecta will move off into space to be detected by our instruments on Earth.Therefore, only the ejecta is considered when comparing to the yields inferred from observations.This is a reasonable assumption, because Arnett's model is only applied to data from the early light curve at maximum light, and it has been demonstrated that surface emission from the remnant only affects the late time evolution of the light curve (Shen & Schwab 2017).
As shown in Figure 14, our simulation produces 56 Ni yields consistent with only the dimmest of the Type Iax population -SN 2019gsc (Srivastav et al. 2020), SN 2008ha (Foley et al. 2009), SN 2010ae (Narayan et al. 2011), SN 2020kyg (Srivastav et al. 2022), SN 2021fcg (Karambelkar et al. 2021), and SN 2007qd (McClelland et al. 2010).These events have ejected 56 Ni masses ranging from 0.0014 -0.013 M ⊙ , while our simulation's range reahes an order of magnitude lower from 0.00030 -0.01528 M ⊙ .Lach et al. (2022)'s C/O WD deflagration study is unable to reach the yields of the very dimmest events, while Kromer et al. (2015)'s C/O/Ne deflagration study comes close (0.0034 M ⊙ ejected 56 Ni).However, as with other studies of deflagrations, the ejected masses produced by our simulation (0.012 -0.133 M ⊙ ) are generally about 10 times lower that that inferred from the dimmest observational data.With the exception of the unusually large range calculated for SN 2021fcg, ejected masses are typically ≥ 0.1 M ⊙ .This relationship is shown in Figure 14. a Ejecta velocity was assumed as this value in order to perform calculations using Arnett's model (Arnett 1982;Valenti et al. 2007).
b These studies perform simulations of deflagrations.For more detail, see Section 1.4.
When comparing the ejected KE and ejected masses of our simulation to observations (Figure 15), it appears that our simulation reproduces the correct ratio, but at too low values.The question is then, is there a way to keep the 56 Ni ejected mass the same, while increasing the total ejected mass and ejected KE of the system in the same ratio?This will be explored in a future 3D study, as discussed further in Section 5.

Core Reignition
It was noticed that in some realizations, during contraction of the central mass and fallback of material onto it, the remnant reached a high enough temperature to reignite.This is why we decided to manually turn off burning at 5s -the burning model described in Section 2.2.2 is not proper for modeling a core reignition.However, we decided to run a few of these realizations with only thermal burning (ϕ CC ) turned on after 3s.There are realizations that detonate, and ones that don't.This brings the question as to how many realizations detonate, and what drives this decision.This is a difficult question to answer generally, and we do not attempt to answer it here, especially given the limits of our 2D simulations.However, this will be very much explored in our future 3D simulation suite, discussed in Section 6.

DISCUSSION AND CONCLUSION
This study performs investigative 2D hydrodynamic simulations of a deflagration in a mixed hybrid C/O/Ne WD to test its viability for explaining the dimmest of the Type Iax events.The simulation is performed in FLASH using a custom 3-stage burning model based on the ADR scheme of Khokhlov (1995) that tracks the burning and energy release.The initial hybrid progenitor is the first created by a stellar evolution code (MESA) rather than a parameterized model, which is evolved from its initial ZAMS state through its accretion phase from a He-donor companion.To obtain a sample, we run 30 "realizations" of the same simulation; the only difference between each realization is the number, shape, and amplitude of perturbations on the flame surface that is used to centrally ignite the star.These realizations were chosen to replicate the initial flames in Augustine et al. ( 2019)'s study of the same progenitor undergoing a DDT.
Similar to other studies of deflagrations (Kromer et al. 2015;Leung & Nomoto 2020;Lach et al. 2022), our simulation's end state shows a dense bound remnant surrounded by sparse ejecta, where the majority of the system's mass is unburned.We obtain ejecta yields of: 0.012 -0.133 M ⊙ total mass, 0.00117 -0.04274 M ⊙ IGE, 0.00030 -0.01528 M ⊙ 56 Ni, and 0.107 -1.660 × 10 49 ergs of KE.The remaining 1.23 -1.35 M ⊙ is either in the dense core left behind, or will fall back onto this core to form a bound remnant.More detail of the yields of our simulation can be found in Table 1.We find that our simulation exhibits a characteristic ejecta velocity of 3500 km/s, which falls within the range of ejecta velocities inferred from observations (Jha 2017).We include only the unbound ejecta when comparing to observational data, however ongoing investigation seeks to determine if and how the trapped elements in the remnant will contribute to luminosity and spectroscopy measurements (Shen & Schwab 2017).
Our study's ejected 56 Ni yields cover the range of the 56 Ni yields of the dimmest of the Type Iax population.However, our total ejecta yields, which are about 10 times larger than the 56 Ni yields, are smaller than the estimations provided by observational studies.This mismatch could stem from both the application of Arnett's analysis to observations, and the setup used in our simulation.
It is important to note that when estimating the 56 Ni yields using Arnett's model, there are assumptions made about a constant opacity and a central 56 Ni distribution (Arnett 1982;Valenti et al. 2007).As Arnett's analysis was created for use of explosions that fully disrupt the star, these assumptions may be less relevant for Type Iax events.However, Arnett's model still provides a valuable comparison to observations, and it is unlikely that it would overestimate ejecta yields by a full order of magnitude, while also having accurate 56 Ni ejecta yields.Therefore, it will be important to generate measures such as light curves and spectra (see Section 6) to make a better comparison with observations.Differences in the progenitor (for example, metallicity), as well as different initial conditions of the model flame (for example, radius) may be able to produce a higher ejecta and KE yield in simulation, while keeping the same low 56 Ni mass that our simulation produces.Changing the initial matchhead radius and geometry (see Section 4.4) was unable to achieve this desired outcome, however other initial conditions of the flame, such as placement, will be varied in a future study.However, other studies of 3D deflagrations have explored this parameter space, the results of which will be summarized here.Lach et al. (2022) found that increasing the central density of their parameterized C/O progenitor resulted in an increase in both IGE and KE yields, while the 56 Ni yield increased to a point, but then decreased with increasing central density.This decrease is because electron capture rates are higher at larger densities, which results in a higher production of more stable isotopes of Ni.Lach et al. (2022) also found that changing the metallicity of the initial ZAMS star did not have a systematic effect on the resulting ejecta composition or energy.Based on these results, it may be the case that a hybrid progenitor with a higher central density would have a higher total IGE and KE yield while keeping the 56 Ni yield the same.
As our progenitor is not a parameterized model, but rather a result of a full stellar evolution calculation, adjusting its central density and metallicity would require re-doing the cal-culations described in Section 2.1 with different assumptions.Physically, the central density of a progenitor is a result of its isolation time -the time it takes for the star to get close enough to its companion to commence accretion.The longer this takes, the more time the progenitor has to cool and even crystallize, which increases its central density at the time of ignition (Lesaffre et al. 2006).For us to create a hybrid progenitor with a higher central density, it would be necessary to increase the distance between the WD and its companion in the MESA simulation described in Section 2.1.However, there is no guarantee that doing this will create the desired hybrid WD, as hybrid WD have a tighter parameter space than traditional C/O ones.This would have to be done with much care.
Deflagrations of C/O/Ne WDs could explain the lowest brightness Type Iax events, but it is highly unlikely that this is the only method that produces a Type Iax event; these hybrid WDs are unusual and have a much lower C-fraction in the core, meaning there is less fuel to burn and turn into 56 Ni.It may be that in the single degenerate paradigm, a range of compositions and densities of the initial progenitor may be able to explain the wide variation of Type Iax events, however the double degenerate paradigm should not be overlooked.Notably, as previously mentioned in Section 1.2, the only observed Type Iax event remnant,"Parker's star", is thought to be created by a C/O + O/Ne WD-WD merger (Gvaramadze et al. 2019;Oskinova et al. 2020;Ritter et al. 2021).It has been suggested that studying more of these remnants may provide clues to determine the progenitor system of the event (Ferrand et al. 2021).
It is also uncertain how and why some of our realizations transition from a deflagration to a detonation upon contraction, incinerating the star and potentially producing a normal SN Ia.As the only difference between each realization is the initial ignition, it is interesting that some realizations detonate, while others do not (as far as we can simulate, ∼ 20s).This shows that the initial flame configuration is exceedingly important in determining the end state of the simulation, and must be carefully studied.Two applicable theories of reignition in a deflagration scenario are the pulsationally assisted gravitationally confined detonation (PGCD) (Jordan et al. 2012), and the pulsating reverse detonation (Bravo & García-Senz 2009).Both mechanisms rely on the pulsation of the bound remnant as it contracts and pulls the surrounding material, which contains burning ashes, back towards it.The key difference between the two is where the detonation ignites: a PGCD will ignite the detonation inside the remnant's dense core where the burning ashes mix with the fuel, while a PRD ignites an inward-moving detonation at the interface of the remnant's core and the fallback ashes.Lach et al. (2022) also tested their models in the gravitationally confined detonation (GCD) scenario, where a deflagration wraps around the star and collides, sparking a detonation.This mechanism produced a much higher ejected 56 Ni content (0.257 -1.057 M ⊙ ), with lightcurves and spectra most similar to 91T-like events, rather than Type Ia or Iax.
If reignition does not occur, the remnant will eventually form a compact object that is enriched with IGEs and C-burning ashes.It has been suggested that these object could be metallic hypervelocity 'rogue stars' that wander outside of galaxies as discussed by Hawkins et al. (2015), under the conditions of an asymmetric initial deflagration that gives the progenitor the kinetic energy required to escape its binary system.In 2D, our symmetry does not allow for an asymmetric deflagration, but it is possible and has been seen in a 3D configuration (Kromer et al. 2015;Lach et al. 2022).A potential candidate rogue star that may have originated from this type of progenitor system is LP 40-365, a metallic rogue star with a velocity of roughly 850 km/s that has been studied for decades.Hermes et al. (2021) determined that based on the WD composition, the star likely had some type of thermonuclear origin, however, this does not explain its velocity.In theory, an asymmetric deflagration could explain the origin of hypervelocity stars like LP 40-365.

FUTURE WORK
This 2D study has helped us identify, explore, and mitigate the challenges of evolving a deflagration in a hybrid progenitor, as well as shown the system to be a viable hypothesis for the dimmest Type Iax events that warrants further investigation.The most critical expansion of this study is to perform 3D simulations and to produce light curves to better compare with observational data.At the same time, exploring the parameter space is necessary to quantify the uncertainty in our results.
The initial flame radius is a free parameter in simulation that corresponds to how much the flame has evolved before the simulation is initialized.It is necessary to include this artificially burned initial region, as it is impossible to directly simulate a realistically sized flame (∼ 1 cm scale) in a full-star (∼ 10 8 cm) simulation.It is important to determine the effect that this initial flame size and placement has on our results.Additionally, it will be important to test whether or not our simulations reignite, creating a delayed detonation and a much brighter event than a Type Iax.The relationship between the initial flame size and perturbations and the likelihood of detonation must be investigated.
Preparing a 3D simulation will involve modeling the pre-existing turbulence in the convective core and running the simulation until the ejecta are in free expansion.The former will be practically accomplished by introducing a turbulence field as in Jackson et al. (2014).The latter requires addressing the issue of potentially unphysical results from our EOS in the difficult low-temperature, low-density regions, discussed in Section 4.2.An important step to solving this is to investigate the low-temperature area immediately in front of the burned plumes, to determine whether they are physically motivated or just a simulation artifact.
To produce light curves, Lagrangian tracer particles will be added to the simulation and post-processed in MESA to determine the elemental abundances for the simulation's end state, following the prescription of Townsley et al. (2019).These abundances will then be separated into those for remnant and ejecta.The ejecta abundances will be mapped to a velocity grid and a radiative transfer calculation will be performed to produce light curves and spectra.This will provide a much better comparison to observations than Arnett's model, although radiative transfer software introduce their own layer of assumptions.
More work needs to be done both in simulation and observational astronomy to understand these dim events.Normally, only the ejecta is used to create light curves to compare to observations, however it is unclear if and how the material in the bound remnant contributes.Some studies indicate that the late-time evolution of the light curve is affected by the remnant (Shen & Schwab 2017).It is also unclear if a bound remnant has been observed in the aftermath of a Type Iax event, and this goal is impeded by the presence of the surviving companion in the binary system.Fortunately, there are active observational campaigns directed towards this goal, such as that of Jha et al. (2021).Additionally, hypervelocity stars such as LP40-365 have been hypothesized to have a Type Iax origin (Hermes et al. 2021).Therefore, it is important to obtain as much detail as possible about the remnant that appears in simulations.Due to the difference in their evolution timescales, it is difficult to simultaneously resolve both the rem-nant and ejecta in simulations.Oftentimes the ejecta are focused on, so that light curves can be created for more direct comparison to observations.While evolving the remnant into a compact object in simulation will take work, it is an effort worth pursuing.Therefore, the goal is to run our simulation until the remnant's pulsation is calmer, and then the remnant abundances will be used as input to a stellar evolution calculation in MESA, to determine the makeup of the compact object that is produced in this process.This will allow a comparison to observations from the literature, and perhaps provide information for future observational campaigns.assumed to be fully ionized, though no checks are made to avoid regions where that is known to be inaccurate, as doing so in regions where pressure ionization is present, rather than the usual Saha ionization, does not have a simple well-defined procedure.More broadly applicable equations of state typically use tabulations in such regions, e.g. the EOS used in the MESA (Jermyn et al. 2023).
In the Helmholtz EOS, the Helmholtz free energy, F , and eight of its partial derivatives are tabulated as a function of density and temperature.P , E int , and S are then computed from this lookup table.This table was computed using the Timmes EOS, which was designed to be accurate at the cost of being computationally intensive (Timmes & Arnett 1999).The Timmes EOS calculated η to 18 significant figures using a Newton-Raphson method, evaluated the Fermi-Dirac integrals exactly in IEEE 64-bit arithmetic, and analytically solved for the derivatives of pressure and thermal energy.To keep this accuracy (IEEE 64-bit precision) but reduce the computation time (by over a factor of 100), these results were used to calculate F and its partial derivatives, which were then tabulated (Timmes & Swesty 2000).The granularity of the table affects the accuracy of the result -we used the same table as that in CASTRO (Almgren et al. 2010), which has values for 271 densities and 101 temperatures in the limits 10 −10 < ρ < 10 11 g/cm 3 and 10 4 < T < 10 11 K1 .Although this table was created with Y e = 1 (pure hydrogen), it can be directly applied to any Y e by making the appropriate adaptations based on Y e .
The choice to tabulate F and its derivatives is dictated by the desire for tight thermodynamic consistency and the use of temperature and density as the desired independent variables, therefore making Helmholtz potential the corresponding fundamental Legendre-transformed potential.P , E int , and S can be easily constructed from F and its derivatives: where the value in the subscript is held constant for the thermodynamic partial derivative.
A biquintic polynomial is used to interpolate the table values in 2D (ρ, T ) space between grid cells, to ensure both that thermodynamic consistency is maintained and that the derivatives of P , ϵ int , and S are continuous across the grid.

Figure 1 .
Figure1.Density (top), temperature (middle), and composition (bottom) of the end state of the hybrid C/O/Ne WD after accreting mass from its companion star.The 29 element reaction network was condensed into 4, so it can be used with our burning model in FLASH.Here, 22 Ne is used as a proxy for neutron-rich material, and the accreted He has been added as 20 Ne as they have the same electron-to-baryon ratio.

Figure 4 .
Figure 4. Plots over time of (a) the total burned mass, (b) unbound mass, (c) unbound IGE mass, and (d) unbound 56 Ni mass.The results of the following six simulations are highlighted: R16 (orange), which produced the most ejected mass; R20 (yellow) and R32 (black), which yielded the most and least ejected 56 Ni mass, respectively; R28 (purple), which by far burned the least mass; and R23 (blue) and R25 (fuchsia), two average runs that produced almost the same ejected 56 Ni mass but exhibit different rise times.Colors are in dark to light ordered by ejected 56 Ni mass.To get a sense of the range and localization of the data, also shown are the IQR for each timestep (purple, shaded), the range of the middle 75% of the data (blue, shaded), and the full range of the data (light blue, shaded).

Figure 5 .
Figure 5. Fraction of material burned to IGE (ϕ qn ) of R28 at 2.0s of simulation time.The curvature of the plumes due to the cylindrical geometry of the setup is apparent.

Figure 6 .
Figure 6.A comparison of R32 (left panels, least ejected 56 Ni mass) and R20 (right panels, most ejected 56 Ni mass) over time.The (a) initial flame, (b) the the fraction of material burned to IGE at 3s, (c) the fraction of material burned to IGE at 20s, and (d) the density at 20s, are shown.

Figure 7 .
Figure 7.A comparison of the half radius (left panels) and double radius (right panels) matchhead sizes for R20 over time.The (a) initial flame, (b) the the fraction of material burned to IGE at 3s, (c) the the fraction of material burned to IGE at 20s, and (d) the density at 20s, are shown.In (b), density contours at 10 2 and 10 4 g/cm 3 are also plotted.

Figure 8 .
Figure 8. Plots over time of (a) the total burned mass, (b) total 56 Ni : IGE Ratio, (c) total IGE mass, and (d) total56 Ni mass for the 75 km (dotted lines) and 300 km (solid lines) radius matchheads.Similarly to Figure4, the results of the following four simulations are highlighted: R20 (yellow) and R32 (black), which yielded the most and least ejected 56 Ni mass for the baseline (150 km) radius matchhead, respectively; and R23 (blue) and R25 (fuchsia), two average runs that produced almost the same ejected 56 Ni mass but exhibit different rise times for the baseline (150 km) radius matchhead.Also shown are the IQR for each timestep (purple, shaded), the range of the middle 75% of the data (blue, shaded), and the full range of the data (light blue, shaded).

Figure 9 .Figure 10 .
Figure 9. Corner plot showing correlations between different final integrated quantities and their distributions.All masses are in units of M ⊙ .Shown are the total burned mass (M Burned Total ), the total 56 Ni mass (M Total 56 Ni ), the total IGE mass (M Total IGE ), the ejected 56 Ni mass (M Ej 56 Ni ), the ejected IGE mass (M Ej IGE ), the total ejected mass (M Ej Total ), the bound 56 Ni mass (M Bound 56 Ni ), the bound IGE mass (M Bound IGE ), and the total bound mass (M Bound Total ).

Figure 11 .
Figure 11.The maximum total kinetic energy of each simulation is plotted against both (top) the final unbound mass, and (bottom) the final total burned mass.The tight correlation between the former indicates that the maximum total KE determines how much mass can escape the gravitational pull of the remnant.The colored points show the same feature highlights as Figure 4.

Figure 12 .
Figure 12.Total unbound kinetic energy (in 10 49 ergs) vs the total unbound mass (in M ⊙ ) of each realization.The colored points show the same feature highlights as Figure 4.

Figure 13 .
Figure 13.The final total 56 Ni yield the final total IGE yield of this study with a deflagration explosion mechanism (grey •) compared to Augustine et al. (2019)'s study with a DDT explosion mechanism (grey ×).

Table 1 .
The range and interquartile range (IQR, middle 50%) of the final yields from our simulation suite.
Note-Bound and unbound material are distinguished based on the value of E binding (defined by Equation IQR away from the 75th and 25th percentiles.Most of the star remains bound, and the majority of the ejecta is unburned fuel.Most of the produced 56 Ni and IGE remains bound.

Table 2 .
The range and interquartile range (IQR, middle 50%) of the final yields from the halved (75 km) matchhead radius simulation suite.

Table 4 .
Ejecta yields calculated from observational data and other simulation studies.Data is ordered by increasing 56 Ni mass while grouping studies of the same event.