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

Catherine Feldman Department of Physics and Astronomy
Stony Brook University
Stony Brook, NY 11794-3800, USA  Institute for Advanced Computational Science
Stony Brook University
Stony Brook, NY 11794-3800, USA  Nathanael Gutierrez School of Computer
Science
Georgia Institute of Technology
Atlanta, GA 30332, USA  Ellis Eisenberg Yale University
New Haven, CT 06520, USA  Donald E. Willcox Center for Computational Sciences
and Engineering
Lawrence Berkeley National Lab
Berkeley, CA 94720, USA  Dean M. Townsley Department of Physics and Astronomy
The University of Alabama
Tuscaloosa, AL 35487-0324, USA  Alan C. Calder Department of Physics and
Astronomy
Stony Brook University
Stony Brook, NY 11794-3800, USA  Institute for Advanced Computational Science
Stony Brook University
Stony Brook, NY 11794-3800, USA

(Accepted September 1, 2023)

###### Abstract

The dimmest and most numerous outlier of the Type Ia supernova population,
Type Iax events, 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
multi-physics 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 matchhead used to start the deflagration. As consistent with earlier
work, these realizations produce the familiar hot dense bound remnant
surrounded by sparse ejecta. Our results indicate the majority of the star
remains unburned ($\sim$ 70%) and bound ($>$ 90%). Our realizations produce
total ejecta yields on the order of $10^{-2}$ \- $10^{-1}$ $M_{\sun}$, ejected
${}^{56}{\rm Ni}$ yields on the order of $10^{-4}$ \- $10^{-2}$ $M_{\sun}$,
and ejecta kinetic energies on the order of $10^{48}$ \- $10^{49}$ ergs.
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 ${}^{56}{\rm Ni}$ yields,
but too-small total yields and kinetic energies. Reignition of the remnant is
also seen in some realizations.

Type Ia supernovae (1728), Supernova remnants (1667), White dwarf stars
(1799), Hydrodynamical simulations (767)

††software: MESA version 8118 (https://mesa.sourceforge.net) (Paxton et al.,
2011, 2013, 2015), FLASH version 4.6.2 (https://flash.rochester.edu) (Fryxell
et al., 2000; Calder et al., 2002) with modifications for C/O and C/O/Ne WDs
(http://astronomy.ua.edu/townsley/code), MESA-to-FLASH hydrostatic profile
builder (https://github.com/Flash-Star/mesautils), (Willcox et al., 2016), yt
(https://yt-project.org) (Turk et al., 2010)

##  1 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_{\text{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}{\rm 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.

###  1.1 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,\text{peak}}<$ -13),
reduced ejecta velocities (2000 $<v_{\rm{ej}}<$ 8000 km/s) and explosion
energies ($\sim 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_{\sun}$, $\sim$ 0.003 – 0.3 $M_{\sun}$ of which is ${}^{56}{\rm 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.

###  1.2 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-to-detonation 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}{\rm 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, 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 single-degenerate 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}{\rm 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.

###  1.3 Hybrid Progenitors

White dwarfs are inactive cores — the end state of low mass ($\lesssim$ 8 - 10
$M_{\sun}$) 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_{\sun}$ form
carbon-oxygen (C/O) WDs, while heavier ones (8 - 10.5 $M_{\sun}$) 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}{\rm C}$, ${}^{16}{\rm O}$, and ${}^{20}{\rm 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 decay of ${}^{13}{\rm
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}{\rm Ni}$
and ejecta yields in DDT simulations due to its lower carbon fraction (Willcox
et al., 2016; Augustine et al., 2019).

###  1.4 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, 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) (LEAFS, C/O progenitor),
Kromer et al. (2015) (LEAFS, unmixed C/O/Ne progenitor), Leung & Nomoto (2020)
(Leung et al. (2015), C/O and mixed C/O/Ne progenitors), and Long et al.
(2014) (FLASH, C/O progenitor). Lach et al. (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.

####  1.4.1 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_{\sun}$ 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}{\rm 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}{\rm Ni}$
in the less dense surrounding material, while those with more ignition points
show the majority of ${}^{56}{\rm 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 ($\alpha$, ${}^{12}{\rm C}$, ${}^{16}{\rm 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 $\times$ 109 g/cm3, central temperature of 6 $\times$ 108 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_{\sun}$, remains bound.
Lach et al. (2022) suggest their findings could account for Type Iax events of
intermediate and high brightness: they report ejected masses between 0.014 and
0.301 $M_{\sun}$and ejected ${}^{56}{\rm Ni}$ masses between 0.0058 - 0.092
$M_{\sun}$. 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^{\rm{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_{\sun}$ C/O WD with a central
density of 3 $\times$ 109 g/cm3, a central temperature of 9 $\times$ $10^{9}$
K, and solar metallicity which is added to the model as a ${}^{22}{\rm 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 $\times$ 109
g/cm3), 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_{\sun}$ is ejected , 0.23 - 0.34
$M_{\sun}$ of which is ${}^{56}{\rm 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.

###  1.5 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_{\sun}$
50/50 C/O core surrounded by a 1.1 $M_{\sun}$, 50/47/3 ${}^{16}{\rm
O}$/${}^{20}{\rm Ne}$/${}^{12}{\rm C}$ shell. A 0.1 $M_{\sun}$ accreted layer
of 50/50 C/O is placed on top of this. The central density is 2.9 $\times
10^{9}$ g/cm3. 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_{\sun}$, 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_{\sun}$ bound remnant with 0.014
$M_{\sun}$ of ejecta. The ejecta is comprised of mostly unburned material,
with 6.4 $\times 10^{-3}$ $M_{\sun}$ of IGE, roughly half of which is
${}^{56}{\rm Ni}$. It has a low ejecta energy of 1.8 $\times 10^{48}$ ergs of
kinetic energy. This simulation gives a smaller ${}^{56}{\rm Ni}$-to-mass
ratio than expected by observations, and also too fast rise and decline of the
light curves. However, its ejected ${}^{56}{\rm 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}{\rm
C}$, 50% ${}^{16}{\rm O}$, and 30% ${}^{20}{\rm Ne}$ by mass, and then adds in
solar metallicity as 2.5% ${}^{22}{\rm Ne}$. It has the same central density
(3 $\times$ 109 g/cm3) as their C/O model. They investigate the effect of
changing the central density of the WD (1.00 - 9.00 $\times$ 109 g/cm3) 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.30 $M_{\sun}$, 0.18 - 0.36
$M_{\sun}$ of which is ${}^{56}{\rm Ni}$. These values are on average slightly
smaller than their C/O WD counterparts.

###  1.6 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.

##  2 Methodology

There are two main pieces needed for this study: a model of a hybrid C/O/Ne WD
progenitor 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.

###  2.1 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, 2013, 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}{\rm
C}$(${}^{12}{\rm C}$,p)${}^{23}{\rm 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}{\rm 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_{\sun}$ hybrid C/O/Ne WD, with a 0.4 $M_{\sun}$ C/O
core with a C fraction of 0.345, enclosed by a 0.69 $M_{\sun}$ O/Ne/Na mantle.

There is a $Y_{e}$ gradient between mantle and core (burned and unburned
material), which in this model is $\Delta Y_{e}$ = -0.00235 (Brooks et al.,
2017). The burned material has a smaller $Y_{e}$ due to the beta decay of
${}^{13}{\rm 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_{\sun}$ 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 $\sim$
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_{\sun}$, with a central density of
$2.2\times 10^{9}$ g/cm3 and a central temperature of $8.2\times 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}{\rm C}$, ${}^{16}{\rm O}$, ${}^{20}{\rm Ne}$, and ${}^{22}{\rm
Ne}$. The prescription for doing this is to keep the ${}^{12}{\rm C}$ mass
fraction ($X_{{}^{12}{\rm C}}$) the same, and then use ${}^{22}{\rm Ne}$ as a
proxy for the neutron-rich nuclides by calculating its mass fraction
($X_{{}^{22}{\rm Ne}}$) using $Y_{e}$. The remaining mass fractions of
${}^{16}{\rm O}$ ($X_{{}^{16}{\rm O}}$) and ${}^{20}{\rm Ne}$ ($X_{{}^{20}{\rm
Ne}}$) were kept in the same ratio as in the MESA profile.

$X_{{}^{12}{\rm C}}^{\rm FLASH}=X_{{}^{12}{\rm C}}^{\rm MESA}$ (1a)
$X_{{}^{22}{\rm Ne}}^{\rm FLASH}=\max(22(\frac{1}{2}-Y_{e}),0.0)$ (1b)
$R=X_{{}^{20}{\rm Ne}}^{\rm MESA}/X_{{}^{16}{\rm O}}^{\rm MESA}\\\ $
$X_{{}^{16}{\rm O}}^{\rm FLASH}=\frac{1-X_{{}^{12}{\rm C}}^{\rm
FLASH}-X_{{}^{22}{\rm Ne}}^{\rm FLASH}}{R+1}$ (1c) $X_{{}^{20}{\rm Ne}}^{\rm
FLASH}=RX_{{}^{16}{\rm O}}^{\rm FLASH}$ (1d)

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.

Figure 1: 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}{\rm Ne}$ is used as a proxy
for neutron-rich material, and the accreted He has been added as ${}^{20}{\rm
Ne}$ as they have the same electron-to-baryon ratio.

Three “layers” can be seen, with transitions between them. The first is the
composition-homogenized core convection zone (1.07 $M_{\sun}$, 1.018$\times
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 $\times 10^{8}$ cm, and contains
ashes from He-burning, where the accreted material has compressed and become
hot enough (1.5-2.5 $\times 10^{8}$ K) to burn He. The remaining $10^{-6}$
$M_{\sun}$ layer that extends to the edge of the star (2.462 $\times 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}{\rm Ne}$ in
the accreted layer — in reality, the accreted layer is mostly ${}^{4}{\rm
He}$, but this is not an element in the burning model. Since ${}^{20}{\rm Ne}$
and ${}^{4}{\rm He}$ have the same electron-to-baryon ratio, this should not
change the equilibrium much.

###  2.2 Simulating a Deflagration

The deflagration simulation was performed using a customized version (Calder
et al., 2007; Townsley et al., 2007, 2009, 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.

####  2.2.1 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 additional 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}<\rho<10^{11}$ g/cm3 and $10^{4}<T<10^{11}$ K (Timmes & Swesty,
2000). However, there is a region in $\rho$ -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}\sim 10^{11}$ ergs, which lies where
$T<\sim 10^{6}$ K and $\rho<\sim 10^{3}$ g/cm3. 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.

####  2.2.2 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 advection-diffusion-
reaction (ADR) scheme of Khokhlov (1995) and Vladimirova et al. (2006). The
addition of ${}^{20}{\rm 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 – $\phi_{fa}$ (fuel to ash),
$\phi_{aq}$ (ash to NSQE), and $\phi_{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\geq\phi_{fa}\geq\phi_{aq}\geq\phi_{qn}$.
$\phi_{fa}$ is evolved using an ADR scheme, and $\phi_{aq}$ and $\phi_{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}{\rm Ne}$ did not
require any substantial changes to to the structure of this model.
${}^{20}{\rm Ne}$ 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}{\rm 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}{\rm C}$-fraction, and ${}^{22}{\rm 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).

####  2.2.3 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 ($\phi_{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):

$r(\theta)=150\text{km}+30\text{km}\sum_{\ell=\ell_{min}}^{\ell_{max}}\sum_{m=-\ell}^{+\ell}A_{\ell}Y^{m}_{\ell}(\theta)$
(2)

The values of $\ell$ are taken from Townsley et al. (2009): $\ell_{max}$ was
chosen so that that the perturbations are resolved, and a relatively large
choice for $\ell_{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:

$r(\theta)=150\text{km}+30\text{km}\sum_{\ell=12}^{16}A_{\ell}Y^{0}_{\ell}(\theta)$
(3)

These $A_{\ell}$ 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_{\ell}$.

####  2.2.4 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 low-density “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\times
10^{14}$ dyne/cm2, 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 ${\sim}4\times 10^{14}$ dyne/cm2 to
${\sim}4\times 10^{15}$ dyne/cm2, 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/cm3.

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/cm3 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 (${\sim}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\times 10^{6}$ K. This high temperature
floor introduces a negligible amount of internal energy into our simulations
at later times ($\sim$ 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.

####  2.2.5 Refinement

We have chosen a fluff density on the order of $10^{-4}$ g/cm3 and the ejecta
reaches densities of the order $10^{-2}$ g/cm3 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\leq\phi_{\rm{fluff}}\leq 1$, where $\phi_{\rm{fluff}}=0$ inside the star,
and $\phi_{\rm{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.

###  2.3 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}{\rm 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”).

##  3 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.

$E_{\rm binding}=E_{\rm grav}+E_{\rm kinetic}$ (4)

Here, $E_{\rm grav}$ is the gravitational potential energy, and $E_{\rm
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}{\rm Ni}$, must be calculated using the progress
variables ($\phi_{fa}$, $\phi_{aq}$, $\phi_{qn}$), neutron excess of the fluid
($Y_{e}$), and initial, unburned mass fractions of ${}^{12}{\rm C}$
($X_{{}^{12}{\rm C}}$), ${}^{20}{\rm Ne}$ ($X_{{}^{20}{\rm Ne}}$), and
${}^{22}{\rm Ne}$ ($X_{{}^{22}{\rm Ne}}$) the proxy for neutron-rich
material).

The simulation considers the remaining material as ${}^{16}{\rm O}$.

$X_{{}^{16}{\rm O}}=1-X_{{}^{12}{\rm C}}-X_{{}^{20}{\rm Ne}}-X_{{}^{22}{\rm
Ne}}$ (5)

The mass of material in each grid cell was calculated using the mass fraction,
density, and 3D volume of each cell.

$M_{i}=X_{i}\rho V$ (6)

where $M_{i}$ is the mass of substance $i$, $X_{i}$ is the mass fraction of
substance $i$, $\rho$ is the density and $V$ the 3D annular volume of a cell.

The mass fraction of IGE is simply the progress variable $\phi_{qn}$, as
described in Section 2.2.2, since the difference between baryon density and
mass density is unsubstantial.

$X_{IGE}=\phi_{qn}$ (7)

The calculation of the mass fraction of ${}^{56}{\rm Ni}$ present in each grid
cell ($X_{{}^{56}{\rm Ni}}$) follows the estimation of Equation 29 of Townsley
et al. (2009), which assumes that ${}^{54}{\rm Fe}$ and ${}^{58}{\rm Ni}$ are
the next most-common elements in the fully burned material.

$X_{{}^{56}{\rm Ni}}=\max\left(\phi_{qn}\frac{Y_{e,n}-Y_{e,{}^{54}{\rm
Fe}+{}^{58}{\rm Ni}}}{Y_{e,{}^{56}{\rm Ni}}-Y_{e,{}^{54}{\rm Fe}+{}^{58}{\rm
Ni}}},0\right)$ (8)

where $Y_{e,i}$ is the electron fraction of substance $i$. The electron
fraction in the fully burned material is

$Y_{e,n}=\frac{Y_{e}-Y_{e,f}(1-\phi_{qn})}{\phi_{qn}}$ (9)

and the electron fraction in the unburned fuel is

$\displaystyle Y_{e,f}=X_{{}^{12}{\rm C}}Y_{e,{}^{12}{\rm C}}+X_{{}^{20}{\rm
Ne}}Y_{e,{}^{20}{\rm Ne}}$ (10) $\displaystyle+X_{{}^{22}{\rm
Ne}}Y_{e,{}^{22}{\rm Ne}}+X_{{}^{16}{\rm O}}Y_{e,{}^{16}{\rm O}}$

and the electron fractions for the remaining substances are calculated as
usual, e.g. $Y_{e,{}^{12}{\rm C}}$ = $Y_{e,{}^{20}{\rm Ne}}$ =
$Y_{e,{}^{16}{\rm O}}$ = 0.5; $Y_{e,{}^{22}{\rm Ne}}=\frac{10}{22}\sim
0.4545$; and $Y_{e,{}^{54}{\rm Fe}+{}^{58}{\rm
Ni}}=\frac{1}{2}(\frac{26}{54}+\frac{28}{58})\sim 0.48212$.

##  4 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).

###  4.1 Overview of Explosion Stages

Figure 2 shows both the density and percentage of material burned to IGE
($\phi_{qn}$) of an average realization, R25, at 0, 0.7, 3, 13, and 20s of
simulation time.

Figure 2: Density (top row) and $\phi_{qn}$ (bottom row) of 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
($\Delta 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\times
10^{7}$ g/cm3 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 $\phi_{aq}$ and $\phi_{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.

###  4.2 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 simulations 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.

###  4.3 Final Yields

Table 1: The range and interquartile range (IQR, middle 50%) of the final
yields from our simulation suite.

Calculated Quantity |  | Total | Bound | Unbound
---|---|---|---|---
Mass ($M_{\sun}$) | Range | 1.36 | 1.230 - 1.351 | 0.012 - 0.133
| Middle 50% | 1.36 | 1.289 - 1.332 | 0.031 - 0.074
Estimated ${}^{56}{\rm Ni}$ ($M_{\sun}$) | Range | 0.066 - 0.127 | 0.064 - 0.115 | 0.00030 - 0.01528
| Middle 50% | 0.095 - 0.110 | 0.090 - 0.100 | 0.00187 - 0.00738
IGE ($M_{\sun}$) | Range | 0.19154 - 0.28051 | 0.18288 - 0.26143 | 0.00117 - 0.04274
| Middle 50% | 0.22135 - 0.25203 | 0.20828 - 0.23183 | 0.00614 - 0.01833
NSQE ($M_{\sun}$) | Range | 0.03286 - 0.05497 | 0.03203 - 0.5497 | 0.00017 - 0.00614
| Middle 50% | 0.04600 - 0.05302 | 0.04454 - 0.04889 | 0.00087 - 0.00298
Ash ($M_{\sun}$) | Range | 0.09079 - 0.16329 | 0.08939 - 0.16149 | 0.00033 - 0.01227
| Middle 50% | 0.11888 - 0.13720 | 0.11610 - 0.13252 | 0.00162 - 0.00647
Fuel ($M_{\sun}$) | Range | 0.90380 - 1.04757 | 0.82241 - 1.01695 | 0.00890 - 0.10458
| Middle 50% | 0.92959 - 0.96152 | 0.87219 - 0.94086 | 0.01783 - 0.05095
${}^{56}{\rm Ni}$ : IGE Ratio | Range | 0.344 - 0.471 | 0.350 - 0.499 | 0.161 - 0.476
| Middle 50% | 0.402 - 0.441 | 0.408 - 0.449 | 0.297 - 0.406
${}^{56}{\rm Ni}$ : Mass Ratio | Range | 0.04848 - 0.09305 | 0.04842 - 0.09274 | 0.02428 - 0.12240
| Middle 50% | 0.0699 - 0.0809 | 0.06880 - 0.08024 | 0.05695 - 0.09506
KE (1049ergs) | Range | 0.89597 - 3.46294 | 0.45480 - 2.29050 | 0.10703 - 1.66041
| Middle 50% | 1.6212 - 2.3807 | 0.9982 - 1.5614 | 0.3445 - 0.8257

Note. — Bound and unbound material are distinguished based on the value of
$E_{\rm binding}$ (defined by Equation 4) having a negative or positive value,
respectively. The IGE, NSQE, Ash, and Fuel yields correspond to
$\mathcal{X}_{N}$, $\mathcal{X}_{q}$, $\mathcal{X}_{a}$, and $\mathcal{X}_{f}$
described in Appendix A.1 and Section 2.2.2, and the estimate of ${}^{56}{\rm
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).

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 at 1.5 $\times$ 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.

Figure 3: A visual representation of selected final integrated quantities from
Table 1. Shown are box-and-whisker plots of (from left to right) the total
ejected mass ($M^{\rm{Ej}}_{\rm{Total}}$), the ejected fuel
($M^{\rm{Ej}}_{\rm{Fuel}}$), the ejected IGE mass ($M^{\rm{Ej}}_{\rm{IGE}}$),
the ejected ${}^{56}{\rm Ni}$ mass ($M^{\rm{Ej}}_{{}^{56}{\rm Ni}}$), the
total ${}^{56}{\rm Ni}$ mass ($M^{\rm{Total}}_{{}^{56}{\rm Ni}}$), the bound
${}^{56}{\rm Ni}$ mass ($M^{\rm{Ej}}_{{}^{56}{\rm Ni}}$)the total IGE mass
($M^{\rm{Total}}_{\rm{IGE}}$), and the bound IGE mass
($M^{\rm{Bound}}_{\rm{IGE}}$). The individual points are data that lie outside
the whisker range, 1.5 $\times$ 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}{\rm Ni}$ and IGE remains bound.

Notably, $\sim$ 66 - 77% of the star remains unburned; $\sim$ 14 - 21% is
fully burned to IGE; only $\sim$ 2 - 4% is partially burned to NSQE; and
$\sim$ 7 - 12% undergoes carbon burning and then stops (the ash state). Only
$\sim$ 5 - 9% of the star is burned to ${}^{56}{\rm Ni}$, which is $\sim$ 40%
of the IGE. In addition, the vast majority, $\sim$ 90 - 99 %, of the star
remains bound, demonstrating that a deflagration does indeed produce much
smaller amounts of ejected material and high-mass elements 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}{\rm 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}{\rm 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}{\rm 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.

Figure 4: Plots over time of (a) the total burned mass, (b) unbound mass, (c)
unbound IGE mass, and (d) unbound ${}^{56}{\rm 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}{\rm 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}{\rm Ni}$ mass but exhibit different
rise times. Colors are in dark to light ordered by ejected ${}^{56}{\rm 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).

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}{\rm 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}{\rm 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.

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

###  4.4 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_{\ell}$.

Figure 6: A comparison of R32 (left panels, least ejected ${}^{56}{\rm Ni}$
mass) and R20 (right panels, most ejected ${}^{56}{\rm 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: 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 102 and 104 g/cm3 are also plotted.

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}{\rm 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}{\rm 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.

Figure 8: Plots over time of (a) the total burned mass, (b) total ${}^{56}{\rm Ni}$ : IGE Ratio, (c) total IGE mass, and (d) total ${}^{56}{\rm Ni}$ mass for the 75 km (dotted lines) and 300 km (solid lines) radius matchheads. Similarly to Figure 4, the results of the following four simulations are highlighted: R20 (yellow) and R32 (black), which yielded the most and least ejected ${}^{56}{\rm 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}{\rm 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). Table 2: The range and interquartile range (IQR, middle 50%) of the final yields from the halved (75 km) matchhead radius simulation suite. Calculated Quantity |  | Total | Bound | Unbound
---|---|---|---|---
Mass ($M_{\sun}$) | Range | 1.36 | 0.993 - 1.296 | 0.069 - 0.372
| Middle 50% | 1.36 | 1.124 - 1.268 | 0.097 - 0.240
Estimated ${}^{56}{\rm Ni}$ ($M_{\sun}$) | Range | 0.113 - 0.159 | 0.094 - 0.126 | 0.007 - 0.040
| Middle 50% | 0.124 - 0.142 | 0.108 - 0.113 | 0.016 - 0.032
IGE ($M_{\sun}$) | Range | 0.239 - 0.314 | 0.190 - 0.246 | 0.016 - 0.079
| Middle 50% | 0.249 - 0.284 | 0.206 - 0.236 | 0.037 - 0.059
${}^{56}{\rm Ni}$ : IGE Ratio | Range | 0.444 - 0.566 | 0.463 - 0.571 | 0.344 - 0.551
| Middle 50% | 0.472 - 0.512 | 0.472 - 0.517 | 0.441 - 0.499
${}^{56}{\rm Ni}$ : Mass Ratio | Range | 0.083 - 0.117 | 0.074 - 0.127 | 0.079 - 0.190
| Middle 50% | 0.091 - 0.104 | 0.089 - 0.099 | 0.105 - 0.157
KE (1049ergs) | Range | 1.901 - 6.041 | 0.867 - 2.119 | 0.711 - 4.666
| Middle 50% | 2.624 - 4.520 | 1.114 - 1.515 | 1.552 - 3.193
Table 3: The range and interquartile range (IQR, middle 50%) of the final
yields from the doubled (300 km) matchhead radius simulation suite.

Calculated Quantity |  | Total | Bound | Unbound
---|---|---|---|---
Mass ($M_{\sun}$) | Range | 1.36 | 1.32815 - 1.34920 | 0.00048 - 0.02142
| Middle 50% | 1.36 | 1.34870 - 1.34895 | 0.00066 - 0.00092
Estimated ${}^{56}{\rm Ni}$ ($M_{\sun}$) | Range | 0.053 - 0.074 | 0.053 - 0.074 | 0.000 - 0.002
| Middle 50% | 0.063 - 0.071 | 0.063 - 0.071 | 0.000 - 0.000
IGE ($M_{\sun}$) | Range | 0.204 - 0.232 | 0.204 - 0.229 | 0.000 - 0.004
| Middle 50% | 0.218 - 0.226 | 0.218 - 0.226 | 0.000 - 0.000
${}^{56}{\rm Ni}$ : IGE Ratio | Range | 0.259 - 0.320 | 0.259 - 0.319 | 0.327 - 0.353
| Middle 50% | 0.288 - 0.314 | 0.288 - 0.314 | 0.333 - 0.346
${}^{56}{\rm Ni}$ : Mass Ratio | Range | 0.039 - 0.054 | 0.039 - 0.055 | 0.000 - 0.071
| Middle 50% | 0.046 - 0.052 | 0.046 - 0.053 | 0.000 - 0.000
KE (1049ergs) | Range | 1.197 - 3.589 | 1.189 - 3.581 | 0.008 - 0.220
| Middle 50% | 2.239 - 3.017 | 2.227 - 2.938 | 0.011 - 0.014

Note. — Eight of the ten realizations produced no ejected mass. The two that
did produced very little, most of which is material from the expanded, unburnt
star.

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
material 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, 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}{\rm Ni}$, and IGE yields decrease, as does the
${}^{56}{\rm 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 of these yields have some overlap with those for the 150 km radius
matchhead (see Table 1) discussed in the rest of this work, with the 75 km
matchhead overlapping the higher range, and the 150 km matchhead overlapping
the lower range. An important observation from this data 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.

###  4.5 Correlations

Figure 9: Corner plot showing correlations between different final integrated
quantities and their distributions. All masses are in units of $M_{\sun}$.
Shown are the total burned mass ($M^{\rm{Burned}}_{\rm{Total}}$), the total
${}^{56}{\rm Ni}$ mass ($M^{\rm{Total}}_{{}^{56}{\rm Ni}}$), the total IGE
mass ($M^{\rm{Total}}_{\rm{IGE}}$), the ejected ${}^{56}{\rm Ni}$ mass
($M^{\rm{Ej}}_{{}^{56}{\rm Ni}}$), the ejected IGE mass
($M^{\rm{Ej}}_{\rm{IGE}}$), the total ejected mass
($M^{\rm{Ej}}_{\rm{Total}}$), the bound ${}^{56}{\rm Ni}$ mass
($M^{\rm{Bound}}_{{}^{56}{\rm Ni}}$), the bound IGE mass
($M^{\rm{Bound}}_{\rm{IGE}}$), and the total bound mass
($M^{\rm{Bound}}_{\rm{Total}}$).

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 includes ${}^{56}{\rm 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}{\rm Ni}$ show a strong positive correlation,
meaning that a reasonably well-defined ${}^{56}{\rm Ni}$:IGE ratio can be
calculated. This ratio is $\sim$ 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_{\sun}$. Unbound quantities tend to
exhibit a stronger correlation than bound quantities, for example the graph of
unbound ${}^{56}{\rm Ni}$ and unbound IGE is more strongly correlated than
that of bound ${}^{56}{\rm 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.

###  4.6 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.

Figure 10: The total KE vs time, with the same feature highlights as Figure 4.
The total KE peaks right as the bu ends, and then decreases as the
gravitational pull of the remnant slows down the star’s expansion. The KE
peaks again once the remnant starts contracting and material falls back on it,
and then decreases again, showing the pulsation of the remnant.

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.

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

Figure 12: Total unbound kinetic energy (in $10^{49}$ ergs) vs the total
unbound mass (in $M_{\sun}$) of each realization. The colored points show the
same feature highlights as Figure 4.

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.

###  4.7 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-by-realization comparison, however a general comparison between
the two explosion mechanisms can certainly be made.

Figure 13: The final total ${}^{56}{\rm Ni}$ yield vs the final total IGE
yield of this study with a deflagration explosion mechanism (grey $\bullet$)
compared to Augustine et al. (2019)’s study with a DDT explosion mechanism
(grey $\times$).

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}{\rm 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}{\rm Ni}$:IGE ratio of the DDT simulation is around 0.75, nearly twice
than that of our deflagration simulation, 0.4.

###  4.8 Comparison to Observations

Using Arnett’s model (Arnett, 1982) or Valenti et al. (2007)’s variation of it
for a inner high-density 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.30 $M_{\sun}$ ejected ${}^{56}{\rm Ni}$.

Table 4: Ejecta yields calculated from observational data and other simulation
studies. Data is ordered by increasing ${}^{56}{\rm Ni}$ mass while grouping
studies of the same event.

Reference | Type Iax | Magnitude | $v_{\rm{ej}}$ | Mass | ${}^{56}{\rm Ni}$ | KE
---|---|---|---|---|---|---
|  |  | km/s | $M_{\sun}$ | $M_{\sun}$ |  $10^{49}$ ergs
Srivastav et al. (2020) | SN 2019gsc | ($M_{r}$) -13.97 $\pm$ 0.16 | 3500 | 0.13 - 0.22 | 0.0014 - 0.0024 | 1 - 2
Tomasella et al. (2020) | SN 2019gsc | ($M_{r}$) -14.28 $\pm$ 0.15 | 3500 | 0.13 | 0.003 $\pm$ 0.001 | 1.2
Foley et al. (2009) | SN 2008ha | ($M_{V}$) -14.21 $\pm$ 0.15 | 2000 | 0.15 | 0.003 $\pm$ 0.0009 | 0.23
Valenti et al. (2009) | SN 2008ha | ($M_{R}$) -14.5 $\pm$ 0.3 | 2300 | 0.1 - 0.5 | 0.003 - 0.005 | 1 - 5
Srivastav et al. (2020) | SN 2008ha |  |  | 0.1 - 0.2 | 0.003 | 1
Stritzinger et al. (2014) | SN 2008ha | ($M_{B}$) -13.79 $\pm$ 0.14 | 4500 - 5500 | 0.4 - 0.5 | 0.004 | 8 - 18
Srivastav et al. (2020) | SN 2010ae |  | 3500a | 0.2 - 0.3 | 0.003 - 0.004 | 3 - 5
Stritzinger et al. (2014) | SN 2010ae | ($M_{B}$) -13.79 $\pm$ 0.14 | 5000 - 6000 | 0.30 - 0.60 | 0.003 - 0.007 | 4 - 30
Srivastav et al. (2022) | SN 2020kyg | ($M_{r}$) -14.91 $\pm$ 0.15 | 3500a | 0.36${}^{+0.08}_{-0.06}$ | 0.007 $\pm$ 0.001 | 4.4${}^{+2.2}_{-1.5}$
Karambelkar et al. (2021) | SN 2021fcg | ($M_{r}$) -12.66 $\pm$ 0.20 |  | 0.05 - 0.4 | 0.008${}^{+0.004}_{-0.005}$ |
McClelland et al. (2010) | SN 2007qd | ($M_{B}$) -15.4 $\pm$ 0.2 | 2800 | 0.15 | 0.013 $\pm$ 0.002 | 2.0
Barna et al. (2020) | SN 2019muj | ($M_{V}$) -16.4 | 5000a | 0.17 | 0.031 $\pm$ 0.005 | 4.4
Magee et al. (2016) | SN 2015H | ($M_{r}$) -17.27 $\pm$ 0.07 | 5000-6000 | 0.5 | 0.06 |
Tomasella et al. (2016) | SN 2014ck | ($M_{B}$) -17.37 $\pm$ 0.15 | 3000 | 0.2-0.5 | 0.10${}^{+0.04}_{-0.03}$ |
Singh et al. (2022) | SN 2020rea |  | 6500a | 0.77${}^{+0.11}_{-0.21}$ | 0.13$\pm$ 0.01 | 19${}^{+2}_{-6}$
Stritzinger et al. (2015) | SN 2005hk | ($M_{B}$) -18.00 $\pm$ 0.25 | 5000 - 7000 | 1.5-2.0 | 0.15-0.20 | 40 - 120
Szalai et al. (2015) | SN 2011ay | ($M_{B}$) -18.15 $\pm$ 0.17 | 6000 | 0.8 | 0.225 $\pm$ 0.010 | 10
Stritzinger et al. (2015) | SN 2012Z | ($M_{B}$) -18.27 $\pm$ 0.09 | 7000 - 9000 | 1.4-2.6 | 0.25-0.29 | 70 - 280
Kromer et al. (2015)b | C/O/Ne WD | ($M_{B}$) -13.2 - -14.6 |  | 0.014 | 0.00340 | 0.18
Lach et al. (2022)b | C/O WD | ($M_{R}$) -15.2 - -17.94 |  | 0.014 - 0.301 | 0.0058 - 0.092 | 0.18 - 9.7
Leung & Nomoto (2020)b | C/O/Ne WD |  | 5000 - 8000 | 0.96 - 1.30 | 0.18 - 0.36 |
Leung & Nomoto (2020)b | C/O WD |  | 5000 - 8000 | 0.92 - 1.36 | 0.23 - 0.34 |

Note. — Empty values indicate this quantity was not calculated.

aEjecta velocity was assumed as this value in order to perform calculations
using Arnett’s model (Arnett, 1982; Valenti et al., 2007).

bThese studies perform simulations of deflagrations. For more detail, see
Section 1.4.

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}{\rm 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}{\rm Ni}$ masses ranging from 0.0014 - 0.013 $M_{\sun}$, while our
simulation’s range reahes an order of magnitude lower from 0.00030 - 0.01528
$M_{\sun}$. 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_{\sun}$ ejected ${}^{56}{\rm Ni}$).

However, as with other studies of deflagrations, the ejected masses produced
by our simulation (0.012 - 0.133 $M_{\sun}$) 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 $\geq$ 0.1 $M_{\sun}$. This relationship is shown in
Figure 14.

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}{\rm 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.

Figure 14: The final ejected mass vs the final ejected ${}^{56}{\rm Ni}$ mass
of this study (grey $\bullet$) compared to observations (blue $\times$) and
other studies: Kromer et al. (2015) (fuscia $\triangledown$) and Lach et al.
(2022) (pink $\triangle$). Figure 15: The final ejected KE vs the final
ejected mass of this study (grey $\bullet$) compared to observations (blue
$\times$) and other studies: Kromer et al. (2015) (fuscia $\triangledown$) and
Lach et al. (2022) (pink $\triangle$).

###  4.9 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 ($\phi_{\rm{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.

##  5 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_{\sun}$ total
mass, 0.00117 - 0.04274 $M_{\sun}$ IGE, 0.00030 - 0.01528 $M_{\sun}$
${}^{56}{\rm Ni}$, and 0.107 - 1.660 $\times$ 1049 ergs of KE. The remaining
1.23 - 1.35 $M_{\sun}$ 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}{\rm Ni}$ yields cover the range of the
${}^{56}{\rm Ni}$ yields of the dimmest of the Type Iax population. However,
our total ejecta yields, which are about 10 times larger than the ${}^{56}{\rm
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}{\rm Ni}$ yields
using Arnett’s model, there are assumptions made about a constant opacity and
a central ${}^{56}{\rm 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}{\rm 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}{\rm 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}{\rm 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}{\rm 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 calculations 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}{\rm
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, $\sim$ 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}{\rm Ni}$ content (0.257 - 1.057 $M_{\sun}$), 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.

##  6 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 ($\sim$ 1 cm
scale) in a full-star ($\sim 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 remnant 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.

The authors thank the anonymous referee, whose input greatly enhanced the
organization of this work. The authors also thank Ken Shen for many rewarding
conversations and suggestions that improved this work.

This work was supported in part by the US Department of Energy under grant DE-
FG02-87ER40317 and by the National Science Foundation under grant 1927880.
Support also came from the Data + Computing = Discovery summer REU program at
the Institute for Advanced Computational Science at Stony Brook University,
supported by the NSF under award 1950052. This work was also supported in part
by the Simons Summer Research Program at Stony Brook University.

The authors would like to thank Stony Brook Research Computing and
Cyberinfrastructure, and the Institute for Advanced Computational Science at
Stony Brook University for access to the high-performance SeaWulf computing
system, which was made possible by a $1.4M National Science Foundation grant
(#1531492).

The software used in this work was developed in part by the DOE NNSA ASC- and
DOE Office of Science ASCR-supported Flash Center for Computational Science at
the University of Chicago.

This research has made use of NASA’s Astrophysics Data System Bibliographic
Services.

##  Appendix A Overview of the Burning Model and Equation of State

###  A.1 Burning Model

In lieu of a full reaction network, the burning is separated into four groups
of baryons ($\mathcal{X}_{i}$) that evolve in three stages ($\phi_{ij}$) as
well as a dynamic final NSE state, the choice, energetics, and timescales of
which are determined by both detailed “self-heating” models – nuclear burning
calculations performed at either constant pressure (isobaric) or constant
density (isochoric) – and calculations of steady-state detonations. These
self-heating calculations used a 200 nuclide reaction network, and included
the effects of electron screening and weak interactions (Calder et al., 2007).
The four groups of baryons are separated by their characteristic reactions: C
and O fuel ($\mathcal{X}_{f}$); C-fusion ash ($\mathcal{X}_{a}$); Si-group
element dominated nuclear statistical quasi-equilibrium (NSQE)
($\mathcal{X}_{q}$); and Fe-group element (IGE) dominated nuclear statistical
equilibrium (NSE) ($\mathcal{X}_{N}$). The transitions between these four
groups are described by corresponding $\phi$ variables: $\phi_{fa}$ tracks
C-burning; $\phi_{aq}$ the conversion of C-burning ash to NSQE; and
$\phi_{qn}$ the relaxation of NSQE to NSE. The properties of the final
material are then evolved to account for changes due to adjustment of the NSE
and electron capture. Here, $\mathcal{X}_{i}$ describe the baryon fraction of
group $i$, and $\phi_{ij}$ describes the fraction of baryons that have
undergone transition from state $i$ to the next consecutive state $j$
($0\leq\phi_{ij}\leq 1$). These are related by Equation 5 of Townsley et al.
(2016):

$\phi_{fa}=\mathcal{X}_{a}+\mathcal{X}_{q}+\mathcal{X}_{N},\phi_{aq}=\mathcal{X}_{q}+\mathcal{X}_{N},\phi_{qn}=\mathcal{X}_{N}$
(A1)

with $1\geq\phi_{fa}\geq\phi_{aq}\geq\phi_{qn}$.

These baryon fractions, rather than mass fractions, are evolved both in the
burning model and in the hydrodynamics solvers, since baryon number and mass-
energy are conserved in fusion reactions, rather than mass. Therefore, when
discussing density, it is really the baryon density that is being referred to:

$\rho\equiv m_{u}n_{B}$ (A2)

where $m_{u}$ is the atomic mass unit and $n_{B}$ is the baryon number
density. However, the baryon density and mass density are similar enough in
value to be used interchangeably when discussing the results of simulations
(Townsley et al., 2009).

Willcox et al. (2016) performed steady-state detonation calculations to
determine the applicability of this three-stage burning model to an initial
composition that includes ${}^{20}{\rm Ne}$, a necessary nuclide in our hybrid
progenitor. This was accomplished by running 225-nuclide steady-state
Zel’doivch, von Neumann, and During (ZND) detonation models with varying mass
fractions of ${}^{20}{\rm Ne}$ (0.01 - 0.45) – the rest of the composition was
chosen to be 0.02 ${}^{22}{\rm Ne}$, 0.48 ${}^{16}{\rm O}$, and the remainder
${}^{12}{\rm C}$. ${}^{20}{\rm Ne}$ was consumed in the first stage of the
burning, $\phi_{fa}$, with the sole effect being slowing the burning, while
the rest of the stages were unaffected. Therefore, the only change required
was to add ${}^{20}{\rm Ne}$ as an element in the first burning stage, to
accurately represent changes in binding energy release and burning
temperature. The three burning stages plus dynamic ash are described below.

The burning model implemented in FLASH describes the consumption of C-fuel by
considering contributions from both deflagration and thermally-activated
detonation. The evolution of $\phi_{fa}$ is described by Equation 25 of
Townsley et al. (2016):

$\frac{D\phi_{fa}}{Dt}={\rm{max}}(0,r_{\rm{RD}})+r_{\rm{CC}}$ (A3)

$r_{\rm{RD}}$ is the contribution from the deflagration front, which is
evolved by an advection-diffusion-reaction (ADR) equation which in general
form may be written as

$\frac{D\phi}{Dt}=\frac{\partial\phi}{\partial
t}+\nabla\cdot(\vec{v}\phi)=\nabla\cdot(\kappa\nabla\phi)+\frac{1}{\tau}R(\phi)$
(A4)

which describes the advection ($\nabla\cdot(\vec{v}\phi)$), diffusion
($\nabla\cdot(\kappa\nabla\phi)$), and reactions ($R(\phi)$) occuring on
timescale $\tau$, of a substance $\phi$ in a fluid moving with velocity
$\vec{v}$ and diffusivity $\kappa$.

The ADR equation for $\phi_{\rm{RD}}$, which represents the model flame, is
presented in Townsley et al. (2009):

$r_{\rm{RD}}=\frac{D\phi_{\rm{RD}}}{Dt}=\kappa\nabla^{2}\phi_{\rm{RD}}+\frac{1}{\tau}R(\phi_{\rm{RD}})$
(A5)

In this form, the choices for the diffusion ($\kappa$) and reaction timescale
coefficient ($\tau$) (described below) are specified to propagate a laminar
flame at a given speed. $R(\phi)$ is described by a “sharpened” Kolmogorov
Petrovski Piskunov (sKPP) equation (Townsley et al., 2007):

$R(\phi)=\frac{1}{4}f(\phi-\epsilon_{0})(1-\phi+\epsilon_{1})$ (A6)

with $\epsilon_{0}$ = $\epsilon_{1}$ = $10^{-3}$ and $f$ = 1.309, constants
that have been chosen to reduce noise in the solution and give the desired
flame width and propagation speed in combination with other parameters defined
here.

Following Townsley et al. (2007), the diffusion coefficient $\kappa$ is
defined as:

$\kappa=\frac{1}{16}sb\Delta x$ (A7)

and the reaction timescale coefficient $\tau$ as:

$\tau=\frac{1}{16}\frac{b\Delta x}{s}$ (A8)

where $b$ = 3.2, $s$ is the flame speed, and $\Delta x$ is the smallest cell
size of the grid.

$r_{\rm{CC}}$ accounts for the contribution of thermal carbon-carbon burning
(ie in detonation), and is not stored as a variable. It is calculated from
Equation 14 of Townsley et al. (2009):

$r_{\rm{CC}}=\rho
X_{{}^{12}C,f}(1-\phi_{fa})^{2}\frac{1}{12}N_{A}\langle\sigma\nu\rangle_{C+C}(T)$
(A9)

where $\rho$ is the local baryon density defined in Equation A2,
$X_{{}^{12}C,f}$ is the ${}^{12}{\rm C}$ fraction in pure fuel, and
$N_{A}\langle\sigma\nu\rangle_{C+C}(T)$ is the temperature-dependent C-C
fusion rate.

$r_{\rm{CC}}$ is introduced in the evolution of $\phi_{fa}$ to allow for a
detonation to consume partially burned material at the artificially thick
flame front. It is therefore suppressed ($r_{\rm{CC}}=0$) where
$\phi_{\rm{RD}}>10^{-6}\land\phi_{fa}-\phi_{\rm{RD}}<0.1$, to prevent
thermally-activating burning in regions where the artificially thick flame is
active, since the cell temperature in such regions is no longer representative
of the actual burning process. As this work focuses on deflagration,
$r_{\rm{CC}}$ is suppressed most of the time.

The other two $\phi$ variables are evolved at tuned timescales (Townsley et
al., 2016).

$\frac{D\phi_{aq}}{Dt}=\frac{\phi_{fa}-\phi_{aq}}{\tau_{\rm{NSQE}}(T)}$ (A10)

$\frac{D\phi_{qn}}{Dt}=\frac{(\phi_{aq}-\phi_{qn})^{2}}{\tau_{\rm{NSE}}(T)}$
(A11)

The inputs to these calculations, i.e. the temperature-dependent timescales
$\tau$ and the flame speed $s$ were determined as follows.

$\tau_{\rm{NSQE}}$ was defined by Calder et al. (2007) using the
aforementioned “self-heating” calculations, as the time taken for the system
evolve from the point where 90% of the ${}^{12}{\rm C}$ was burned, to the
point of maximum ${}^{28}{\rm Si}$ abundance:

$\tau_{\rm{NSQE}}(T)={\rm{exp}}(\frac{182.06}{T_{9}}-46.054)$ (A12)

$\tau_{\rm{NSE}}$ was tuned to the Si-group consumption timescale of a
200-nuclide, plane-parallel, steady-state Zel’doivch, von Neumann and During
(ZND) detonation model (Townsley et al., 2016, and references therein). This
was accomplished by fitting the curves of Si-group consumption over time
evaluated at pre-shock densities between 0.3 and 10 $\times$ 107 g/cm3.

$\tau_{\rm{NSE}}(T)={\rm{exp}}(\frac{201.0}{T_{9}}-46.77)$ (A13)

The flame speed, $s$, is dependent on the local composition and density. It is
taken to be the maximum of the laminar flame speed, $s_{\rm{lam}}$, and a
minimum flame speed, $s_{\rm{min}}$ that accounts for an increase in flame
speed from the increased surface area of the flame as it undergoes the RT
instability.

$s=\rm{max}(s_{\rm{min}},s_{\rm{lam}})$ (A14)

where

$s_{\rm{min}}=0.5\sqrt{Agm\Delta x}$ (A15)

where $A$ is the C-burning Atwood number (between the unburned fuel and
C-burning ash), $g$ is the acceleration due to gravity, $m$ is a tuneable
parameter, $0.04\leq m\leq 0.06$, and $\Delta x$ is the grid resolution.

$A$ is linearly interpolated from a table generated at the beginning of the
simulation for a range of densities and ${}^{12}{\rm C}$-fractions, using the
global ${}^{22}{\rm Ne}$-fraction. $s_{\rm{lam}}$ is similarly interpolated
from a table of density, ${}^{12}{\rm C}$-fraction, and ${}^{22}{\rm
Ne}$-fraction, which was obtained from detailed reaction network calculations
by Chamulak et al. (2007), augmented by an interpolation of those performed by
Timmes & Woosley (1992) at lower densities (Townsley et al., 2009). For the
hybrid model that includes ${}^{20}{\rm Ne}$, $s_{\rm{lam}}$ is simply
extrapolated to lower densities and carbon fractions. Since the flame
propagation is dominated by the RT instability and turbulence, this is assumed
to be a reasonable approximation (Willcox et al., 2016).

Note for clarity that the above prescription for the flame speed is for 2D
simulations and does not include corrections for turbulence. Future 3D
simulations will use a background turbulent velocity field and appropriate
corrections to speed from the turbulence-flame-interaction (Jackson et al.,
2014).

As seen in the direct numerical simulations of Calder et al. (2007), the
realistic flame is on the order of $10^{-3}$ cm at the densities of the
interior of the star ($10^{9}$g/cm3). It is not computationally possible to
reach this resolution while simultaneously resolving the WD, so the above
description of the flame is practically implemented by artificially thickening
it, so that the full range of $\phi_{fa}$ ([0,1]) always stretches across 4
grid zones. Townsley et al. (2009) showed that the minimum resolution
necessary to accurately describe the model flame was $\Delta x=$4 km. In
reality, as the density of material burned by a flame decreases, the flame
front widens, and this is described in our model by the energy release
timescales $\tau_{\rm{NSQE}}$ and $\tau_{\rm{NSE}}$.

The final piece of the burning model is the adjustment of the properties of
the material in NSQE and NSE, which is referred to as the dynamic ash state.
While corresponding fusion and photo-dissociation (strong reactions) are held
in balance, changes in density and temperature from the bulk fluid motion as
well as neutronization from electron-capture (weak interactions) affect the
equilibrium of this state. This dynamic ash state must therefore capture
changes that occur on both the same timescales as the hydrodynamics and on
longer timescales where weak interactions become important. To achieve this,
the simulation evolves additional variables that track the electron fraction
$Y_{e}$, the average number of nucleons per nucleus $\bar{A}$, and the average
nuclear binding energy per baryon $\bar{q}$, which evolve both during the
burning and in the final NSE state where they become especially important in
influencing the amount of ${}^{56}{\rm Ni}$ produced (Calder et al., 2007).

To determine how these variables change during the dynamic ash state,
Seitenzahl et al. (2009) expanded upon the work of Calder et al. (2007) and
performed calculations with 443 nuclides, which were then turned into tables.
Tables for both isobaric and isochoric calculations are available for 0.44
$\leq Y_{e}\leq$ 0.50125, and using the set of values {$Y_{e}$, $\rho$, $T$},
$\dot{Y}_{e}$ and and the energy loss rate due to neutrinos ($\epsilon_{\nu}$)
can be found by interpolation. From these, the resulting change in $Y_{e}$ and
the rate of change in nuclear binding energy ($\epsilon_{\rm{nuc}}$) of the
NSE state can be calculated, as described in Townsley et al. (2007):

$\dot{Y}_{e,\rm{NSE}}=\phi_{qn}\dot{Y}_{e}(Y_{e},\rho,T)$ (A16)

$\epsilon_{\rm{nuc}}=\frac{\Delta\bar{q}}{\Delta
t}-\phi_{qn}[\dot{Y}_{e}N_{A}c^{2}(m_{p}+m_{e}-m_{n})+\epsilon_{\nu}]$ (A17)

where $N_{A}$ is Avogadro’s number and $m_{p}$, $m_{e}$, and $m_{n}$ are the
mass of the proton, electron, and neutron, respectively.

###  A.2 Equation of State

The EOS used for our application is the Helmholtz EOS described in Timmes &
Swesty (2000) with additional Coulomb corrections. It calculates pressure
($P$), internal energy ($E_{\rm{int}}$), and entropy ($S$) contributions
separately for radiation, ions, electrons, positrons, and Coulomb effect
corrections, and sums them together. For radiation, a thermal photon
distribution is used whose temperature matches that of the other components;
the ions are assumed to act as an ideal gas; the electrons and positrons are
modeled as a non-interacting Fermi-gas; and the Coulomb corrections account
for the interaction between the ions and free electrons. Calculating the state
variables of the Fermi-gas requires evaluating the chemical potential ($\eta$)
using the Fermi-Dirac integrals. The material is 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_{\rm{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 $\eta$ 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}<\rho<10^{11}$ g/cm3 and
$10^{4}<T<10^{11}$ K111This extended table is available at:
https://github.com/AMReX-
Astro/Microphysics/blob/main/EOS/helmholtz/helm_table.dat. 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_{\rm{int}}$,
and $S$ can be easily constructed from $F$ and its derivatives:

$P=\rho^{2}\left(\frac{\partial F}{\partial\rho}\right)_{T}$ (A18)

$E_{\rm{int}}=F+TS$ (A19)

$S=-\left(\frac{\partial F}{\partial T}\right)_{\rho}$ (A20)

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 ($\rho$,
$T$) space between grid cells, to ensure both that thermodynamic consistency
is maintained and that the derivatives of $P$, $\epsilon_{\rm{int}}$, and $S$
are continuous across the grid.

## References

  * Abbott et al. (2019) Abbott, T. M. C., Allam, S., Andersen, P., et al. 2019, 872, L30, doi: 10.3847/2041-8213/ab04fa
  * Almgren et al. (2010) Almgren, A. S., Beckner, V. E., Bell, J. B., et al. 2010, The Astrophysical Journal, 715, 1221, doi: 10.1088/0004-637x/715/2/1221
  * Arnett (1982) Arnett, D. 1982, The Astrophysical Journal, 253, 785, doi: 10.1086/159681
  * Augustine et al. (2019) Augustine, C. N., Willcox, D. E., Brooks, J., Townsley, D. M., & Calder, A. C. 2019, The Astrophysical Journal, 887, 188, doi: 10.3847/1538-4357/ab511a
  * Barna et al. (2020) Barna, B., Szalai, T., Jha, S. W., et al. 2020, Monthly Notices of the Royal Astronomical Society, 501, 1078, doi: 10.1093/mnras/staa3543
  * Bravo & García-Senz (2009) Bravo, E., & García-Senz, D. 2009, The Astrophysical Journal, 695, 1244, doi: 10.1088/0004-637X/695/2/1244
  * Brooks et al. (2016) Brooks, J., Bildsten, L., Schwab, J., & Paxton, B. 2016, The Astrophysical Journal, 821, 28, doi: 10.3847/0004-637x/821/1/28
  * Brooks et al. (2017) Brooks, J., Schwab, J., Bildsten, L., Quataert, E., & Paxton, B. 2017, The Astrophysical Journal, 834, L9, doi: 10.3847/2041-8213/834/2/l9
  * Calder et al. (2010) Calder, A. C., Krueger, B. K., Jackson, A. P., et al. 2010, arXiv e-prints, arXiv:1011.4022.  https://arxiv.org/abs/1011.4022
  * Calder et al. (2002) Calder, A. C., Fryxell, B., Plewa, T., et al. 2002, The Astrophysical Journal Supplement Series, 143, 201, doi: 10.1086/342267
  * Calder et al. (2007) Calder, A. C., Townsley, D. M., Seitenzahl, I. R., et al. 2007, The Astrophysical Journal, 656, 313, doi: 10.1086/510709
  * Chamulak et al. (2007) Chamulak, D. A., Brown, E. F., & Timmes, F. X. 2007, ApJ, 655, L93, doi: 10.1086/511856
  * Colella & Glaz (1985) Colella, P., & Glaz, H. M. 1985, Journal of Computational Physics, 59, 264, doi: 10.1016/0021-9991(85)90146-9
  * De Gerónimo et al. (2022) De Gerónimo, F. C., Miller Bertolami, M. M., Plaza, F., & Catelan, M. 2022, A&A, 659, A150, doi: 10.1051/0004-6361/202142341
  * Denissenkov et al. (2013) Denissenkov, P. A., Herwig, F., Truran, J. W., & Paxton, B. 2013, The Astrophysical Journal, 772, 37, doi: 10.1088/0004-637x/772/1/37
  * Farmer et al. (2015) Farmer, R., Fields, C. E., & Timmes, F. X. 2015, The Astrophysical Journal, 807, 184, doi: 10.1088/0004-637x/807/2/184
  * Ferrand et al. (2021) Ferrand, G., Warren, D. C., Ono, M., et al. 2021, The Astrophysical Journal, 906, 93, doi: 10.3847/1538-4357/abc951
  * Foley et al. (2009) Foley, R. J., Chornock, R., Filippenko, A. V., et al. 2009, The Astronomical Journal, 138, 376, doi: 10.1088/0004-6256/138/2/376
  * Foley et al. (2013) Foley, R. J., Challis, P. J., Chornock, R., et al. 2013, The Astrophysical Journal, 767, 57, doi: 10.1088/0004-637x/767/1/57
  * Fryxell et al. (2000) Fryxell, B., Olson, K., Ricker, P., et al. 2000, The Astrophysical Journal Supplement Series, 131, 273, doi: 10.1086/317361
  * Gvaramadze et al. (2019) Gvaramadze, V. V., Gräfener, G., Langer, N., et al. 2019, Nature, 569, 684, doi: 10.1038/s41586-019-1216-1
  * Hansen et al. (2004) Hansen, C. J., Kawaler, S. D., & Trimble, V. 2004, Stellar Interiors: Physical Principles, Structure, and Evolution, 2nd edn. (Springer Science+Business Media New York)
  * Hawkins et al. (2015) Hawkins, K., Kordopatis, G., Gilmore, G., et al. 2015, Monthly Notices of the Royal Astronomical Society, 447, 2046, doi: 10.1093/mnras/stu2574
  * Hermes et al. (2021) Hermes, J. J., Putterman, O., Hollands, M. A., et al. 2021, The Astrophysical Journal Letters, 914, L3, doi: 10.3847/2041-8213/ac00a8
  * Jackson et al. (2014) Jackson, A. P., Townsley, D. M., & Calder, A. C. 2014, The Astrophysical Journal, 784, 174, doi: 10.1088/0004-637x/784/2/174
  * Jermyn et al. (2021) Jermyn, A. S., Schwab, J., Bauer, E., Timmes, F. X., & Potekhin, A. Y. 2021, The Astrophysical Journal, 913, 72, doi: 10.3847/1538-4357/abf48e
  * Jermyn et al. (2023) Jermyn, A. S., Bauer, E. B., Schwab, J., et al. 2023, The Astrophysical Journal Supplement Series, 265, 15, doi: 10.3847/1538-4365/acae8d
  * Jha (2017) Jha, S. W. 2017, in Handbook of Supernovae, ed. A. W. Alsabti & P. Murdin, 375, doi: 10.1007/978-3-319-21846-5_42
  * Jha et al. (2021) Jha, S. W., Barna, B., Camacho-Neves, Y., et al. 2021, Radioactive Stars: Bound Remnants from White Dwarf Supernovae, HST Proposal. Cycle 29, ID. #16683
  * Jordan et al. (2012) Jordan, G. C., Graziani, C., Fisher, R. T., et al. 2012, The Astrophysical Journal, 759, 53, doi: 10.1088/0004-637X/759/1/53
  * Karambelkar et al. (2021) Karambelkar, V. R., Kasliwal, M. M., Maguire, K., et al. 2021, The Astrophysical Journal Letters, 921, L6, doi: 10.3847/2041-8213/ac2e90
  * Khokhlov (1995) Khokhlov, A. M. 1995, ApJ, 449, 695, doi: 10.1086/176091
  * Kromer et al. (2015) Kromer, M., Ohlmann, S. T., Pakmor, R., et al. 2015, Monthly Notices of the Royal Astronomical Society, 450, 3045, doi: 10.1093/mnras/stv886
  * Lach et al. (2022) Lach, F., Callan, F., Bubeck, D., et al. 2022, A&A, 658, 27, doi: 10.1051/0004-6361/202141453
  * Lach et al. (2022) Lach, F., Callan, F. P., Sim, S. A., & Röpke, F. K. 2022, A&A, 659, A27, doi: 10.1051/0004-6361/202142194
  * Lesaffre et al. (2006) Lesaffre, P., Han, Z., Tout, C. A., Podsiadlowski, P., & Martin, R. G. 2006, MNRAS, 368, 187, doi: 10.1111/j.1365-2966.2006.10068.x
  * Leung et al. (2015) Leung, S.-C., Chu, M.-C., & Lin, L.-M. 2015, Monthly Notices of the Royal Astronomical Society, 454, 1238, doi: 10.1093/mnras/stv1923
  * Leung & Nomoto (2020) Leung, S.-C., & Nomoto, K. 2020, ApJ, 900, 54, doi: 10.3847/1538-4357/aba1e3