# Physics of Thermonuclear Explosions: Magnetic Field Effects on Deflagration
Fronts and Observable Consequences

Boyan Hristov Center for Space Plasma and Aeronomic Research, University of
Alabama in Huntsville, Huntsville, Alabama  Peter Hoeflich Department of
Physics, Florida State University, Tallahassee, Florida  David C. Collins
Department of Physics, Florida State University, Tallahassee, Florida

(Received today; Revised today; Accepted today)

###### Abstract

We present a study of the influence of magnetic field strength and morphology
in Type Ia Supernovae and their late-time light curves and spectra. In order
to both capture self-consistent magnetic field topologies as well evolve our
models to late times, a two stage approach is taken. We study the early
deflagration phase ($\sim 1\rm{s}$) using a variety of magnetic field
strengths, and find that the topology of the field is set by the burning,
independent of the initial strength. We study late time ($\sim 1000$ days)
light curves and spectra with a variety of magnetic field topologies, and
infer magnetic field strengths from observed supernovae. Lower limits are
found to be $10^{6}\rm{G}$. This is determined by the escape, or lack thereof,
of positrons that are tied to the magnetic field. The first stage employs 3d
MHD and a local burning approximation, and uses the code Enzo. The second
stage employs a hybrid approach, with 3D radiation and positron transport, and
spherical hydrodynamics. The second stage uses the code HYDRA. In our models,
magnetic field amplification remains small during the early deflagration
phase. Late-time spectra bear the imprint of both magnetic field strength and
morphology. Implications for alternative explosion scenarios are discussed.

supernovae:general, instabilities, magnetic fields, magneto-hydrodynamics,
turbulence, radiative transfer

††journal: ApJ

##  1 Introduction

Thermonuclear SNe, or  Type Ia (SNe Ia), are explosions of white dwarfs (Hoyle
& Fowler, 1960). They are important for understanding the Universe and the
origin of elements, and are a powerful tool for measuring large distances.
They are also laboratories for understanding the physics of flames,
hydrodynamic instabilities, radiation transport, non-equilibrium systems, and
nuclear and high energy physics in regimes not accessible by ground-based
experiments. Here, we examine the impact of magnetic fields on SNe Ia, which
may alter the sphericity, the nuclear burning front, light curves and spectral
properties.

While SNe Ia light curves can be used as ’quasi-standard candles‘ (Phillips,
1993), there is growing observational evidence for spectral diversity among
SNe Ia that may impact their accuracy as distance measures and, thus, the use
of SNe Ia for precision cosmology. This has prompted sub-classifications based
on observational characteristics, e.g. high- and low-velocity SNe Ia (Benetti
et al., 2005), shallow, core-normal, broad Si-lines, or cool SNe Ia (Branch et
al., 2005). These classifications are widely used for modern data sets (Branch
et al., 2009; Wang et al., 2013; Folatelli et al., 2013). The source of these
spectral differences may come from similar but aspherical objects that are
seen from different angles (Höflich et al., 2006; Motohara et al., 2006; Maeda
et al., 2010; Shen et al., 2018), or they may indicate a differences in
progenitor or explosion scenarios (Hoeflich & Khokhlov, 1996; Quimby et al.,
2006; Shen et al., 2010; Polin et al., 2019). Likely, it is a combination of
both.

A white dwarf, left alone, will eventually cool to the background temperature
of the universe. For the WD to explode, it must interact with a close
companion during the progenitor evolution leading to the explosion. The
companion may either be another WD, (a Double Degenerate System, DD) with
short orbital period, or a non-degenerate star (a Single Degenerate System,
SD) such as a main-sequence, a helium, or a Red-Giant star.

While the exact scenario or scenarios leading to the explosion are still under
study, they can be classified by three distinct mechanisms that trigger the
explosion. These are; slow compressional heating, surface helium detonations,
and the merger of two white dwarfs. The first happens on long timescales,
while the latter two occur quickly. We will discuss each of these in turn.

In the first explosion scenario, the WD accretes material from a companion in
either a DD system on long time scales, so called secular mergers, or in a SD
system (Whelan & Iben, 1973; Piersanti et al., 2003). The explosion is
triggered by compressional heat close to the center of the WD when approaching
the critical Chandrasekhar mass, $M_{Ch}$. The flame propagates by
deflagration (Nomoto K., 1984) and, more likely, starts as a (subsonic)
deflagration and transitions to a (supersonic) detonation ( deflagration-to-
detonation transition (DDT)). The transition from deflagration to detonation
is likely due the mixing of burned and unburned matter, called the Zeldovich
mechanism. (Khokhlov, 1995a; Niemeyer et al., 1996). The mixing process can be
understood in terms of the Zeldovich gradient mechanism (Brooker et al.,
2021); or a unified turbulence-induced mechanism that makes DDT unavoidable
(Poludnenko et al., 2019) at densities suggested by observations. See (Höflich
et al., 2013) for a detailed discussion.

In the second explosion scenario, a surface helium detonation (HeD)111These
have also been referred to as DD for _Double Detonations_ in the literature
triggers a detonation of a sub-$M_{\rm Ch}$WD with a C/O core. (Woosley et
al., 1980; Nomoto, 1982; Livne, 1990; Woosley & Weaver, 1994; Hoeflich &
Khokhlov, 1996; Kromer et al., 2010; Sim et al., 2010; Woosley & Kasen, 2011;
Shen & Moore, 2014; Tanikawa et al., 2018; Glasner et al., 2018; Townsley et
al., 2019). For triggering the initial detonation, these models require an
unmixed surface He layer of $\approx 10^{-2...-1}M_{\odot}$ or $\approx
5\times 10^{-3...-2}M_{\odot}$ with assumed mixing of carbon into the helium
layer on microscopic scales.

In the third possible scenario, two WDs merge or collide. This process occurs
on a dynamical timescale, much faster than the slow accretion timescales from
the previous processes (Webbink, 1984; Iben & Tutukov, 1984; Benz et al.,
1990; Rasio & Shapiro, 1994; Hoeflich & Khokhlov, 1996; Segretain et al.,
1997; Yoon et al., 2007; Wang et al., 2009b, a; Lorén-Aguilar et al., 2009;
Pakmor et al., 2011a; Isern et al., 2011; Pakmor et al., 2012a). In
simulations of this process, the ejecta show large-scale density asymmetries.

It should also be mentioned that any of these triggering mechanisms can be
realized within the common envelope of an Asymptotic Giant Branch (AGB) star
and a WD. This has been invoked to explain super luminous SNeIa \- the so
called Super-$M_{\rm Ch}$explosions. Here, the degenerate C/O core of the AGB
star can explode in two ways. Either, the core grows by accretion on secular
time scales from a disrupted WD and compressional heat triggers the
thermonuclear runaway, or the degenerate AGB core merges with the WD on a
dynamical timescale (Hoeflich & Khokhlov, 1996; Yoon et al., 2007; Kashi &
Soker, 2011; Hoeflich, 2017; Hsiao et al., 2020).

The predominant mechanism for ’normal’ SNe Ia is still under debate with
favorites changing with time. As we will discuss in the conclusions, nuclear
physics dominates the final outcome of the explosion, and successful models
result in almost overlapping progenitor mass ranges among different scenarios.

In our study we examine the impact of magnetic fields on several aspects of
SNe Ia in the classical delayed-detonation scenario because this allows us to
reproduce light curves and spectra of classical SNe Ia. The results will be
put into context of other explosion scenarios in the final discussion and
conclusions.

Our study is based on and combines two previous papers. Hristov et al. (2018)
(Paper I) showed the magnetic field effects on nuclear burning fronts in an
rectangular tube. Penney & Hoeflich (2014) (Paper II) studied the effect of
magnetic fields for dipole and turbulent morphologies on positron transport,
which showed the presence of high magnetic fields ($B$-fields) in some
observed supernovae. The goal of this paper is to extend and combine these
studies for full star simulations with initial conditions which allow to
reproduce observations, and to test the basic assumption on the field
morphology in Paper II.

The obvious question concerns the origin of high $B$-fields. While some WDs
are observed with $B$ of several $10^{7}\,\,\mathrm{G}$, the majority have no
measurable field (Liebert et al., 2003; Schmidt et al., 2003; Silvestri et
al., 2007; Tout et al., 2008). As will be discussed in the conclusion, high
$B$ would require field amplification: a) by rotation induced circulation b)
during the smoldering phase, which is characterized by convection driven non-
explosive carbon burning, with large Eddie sizes corresponding to the pressure
scale height in the WD (Hoeflich & Stein, 2002); or c) during the dynamical
phase of the explosion considered here.

For the DDT mechanism, one of the crucial problems is how to partially
suppress the strong Rayleigh-Taylor (RT) instabilities during the explosion
(Khokhlov, 1995b; Niemeyer & Kerstein, 1997; Gamezo et al., 2003a; Röpke,
2005; Hoeflich, 2006; Motohara et al., 2006; Penney & Hoeflich, 2014; Diamond
et al., 2015; Galbany et al., 2019; Yang et al., 2019). This instability forms
when low-density, high-temperature material is accelerated into higher density
material. The acceleration is most commonly a gravitational field, but can be
also bulk acceleration. The instability manifests as a pattern of rising
mushroom-like structures of low density burned material, leading to large
scale mixing. In particular, 3D hydrodynamical simulations predict rising
plumes throughout the entire WD (Gamezo et al., 2003a; Röpke, 2005). Although
plumes at the RT scales have been observed in supernovae remnants (Fesen et
al., 2007) and indicated by high-resolution polarization spectra of SNe Ia
(Patat et al., 2012; Yang et al., 2019), they are very constrained in velocity
space. High magnetic fields are known to suppress the RT instability
(Chandrasekhar, 1961), and as a WD is a fully ionized plasma, it is reasonable
to consider the effect of magnetic fields (Remming & Khokhlov, 2014; Hristov
et al., 2018).

In Paper I, we simulated a $240\rm{km}\times 15\rm{km}\times 15\rm{km}$ tube
of constant density with WD conditions with magnetic fields from
$10^{9}\rm{G}$ to $10^{12}\rm{G}$. Our finding was that the development of RT
instability was partially suppressed at magnetic field strength of
$10^{10}\rm{G}$, and almost completely suppressed at $10^{12}\rm{G}$. In most
cases, this was related to a reduction in the overall burning rate. However,
certain configurations of burning front and magnetic field (namely models Z12
and YZ12 where the field is strong and aligned with the propagation direction)
the front speed and thickness of the burning region are increased relative to
the unmagnetized case. In this configuration, the field maintains a thicker
burning region, and the increased burning increases the speed of the front. If
this configuration manifests in a real star, it may give another route to the
DDT.

Paper II studied the effect of magnetic fields for dipole and turbulent
morphologies on the positron transport, and found that fields can greatly
impact the late-time light curves and line profiles. The size of the effect
greatly depends on the morphology of the field. In previous studies, the
magnetic fields have been assumed (a) to be radial or locally-trapping small-
scale fields (Milne et al., 2001); or (b) to be dipoles or arbitrary turbulent
fields (Höflich et al., 2004).

The trapping of positrons by magnetic fields can be observed as an increase in
brightness in the late-time light curve, and in the width and evolution of
certain atomic line profiles. Between 300 and 1000 days after the explosion,
the primary energy source in the supernova is positrons produced by the decay
of ${}^{56}{\rm{Co}}$. If the positrons are trapped by the magnetic field, the
supernova will remain brighter than if they are allowed to freely diffuse away
from their birthplace. The degree of trapping depends on both the morphology
and strength of the magnetic field, and has been subject to several PhD theses
(P.A. Milne, 1999; Penny et al., 2012). For light curve studies, Milne et al.
(2001) calculated positron transport for radial magnetic fields and assumed
that, for turbulent fields, all positrons are locally trapped.

In Paper II, we used small scale turbulent and dipole fields as representation
of large fields. We showed that positrons are not locally trapped even for
turbulent fields. The results depend on the size, scale and morphology of the
field. We also showed that magnetic fields and subsequent positron trapping
can also impact late-time near infrared (NIR) spectra. Specifically we focused
on the 1.644 $\mu m$ feature, which is dominated by a single [FeII]
transition, rather than consisting of other blended features in the optical
and NIR (Hoeflich et al., 2004; Diamond et al., 2015). After about day 300,
the shape of the line profile provides important information about the
distribution of ${}^{56}{\rm{Ni}}$, including asymmetries, which may in part
come from magnetic fields (Hoeflich et al., 2004; Motohara et al., 2006;
Diamond et al., 2015; Maguire et al., 2018). Using this line, lower limits for
the magnetic field in SN2003du field were found to be $5000G$ using analytical
approximations for the positron transport (Hoeflich et al., 2004). Detailed
positron transport calculations and observations of NIR line profiles suggest
lower limits larger by one to two orders of magnitude.

The 1.644 $\mu\rm{m}$ line also opened a new aspect and further motivates our
focus on DDT models. While the line offset may be influenced by peculiar
velocity of the progenitor system relative to the host galaxy (typically 100
$\,\,\mathrm{km/s}$ ), it also reflects the orbital velocity of the binary
system. Narrow binary systems result in large offsets in velocity, whereas
wide binary systems (and lower velocity offsets) are expected when the
companion star is non-degenerate. Though the statistics is small, the offset
shift ($v<500...1,000km/sec)$ of the feature in several SNe Ia ( SN2003du
(Hoeflich et al., 2004), SN2005df (Diamond et al., 2014), SN2014J (Diamond et
al., 2018), SN2012ht, SN2013aa, and SN2012cg (Maguire et al., 2018)) suggests
the existence of wide progenitor systems rather than narrow systems to be
expected for single-degenerate progenitors. Two possible exception are
SN2012fr and SN2013cf (Maguire et al., 2018), which show velocities consistent
with double-degenerate progenitors (Shen et al., 2018). Both narrow and wide
systems are can lead to $M_{\rm Ch}$ explosions, whereas mergers and helium
detonation models can only originate from close binaries with large
velocities.

This work extends and connects two prior studies with a focus on the
morphology and size of the magnetic field. Namely, it tries to answer the
following questions: Do the conclusions from simulations in the box still hold
up for full-stars with density gradients? Does a deflagration phase change the
morphology of an initial dipole field? Can the high magnetic fields be
produced during the early deflagration phase? What constraints can be put on
the size of the magnetic fields in the progenitor without arbitrarily assuming
its morphology?

The goal of this paper is to study the effect of magnetic fields on the
explosion and the observational consequences. For computational feasibility,
we do not present a single end-to end simulations for the entire evolution but
strive to perform two interconnected and consistent sets of simulations. In
the first part, we assume the size of the magnetic field and determine its
morphology from simulation. In the second part, we assume the morphology of
the magnetic field and constrain its size from light curves. It is important
to address the question whether observations can constrain the magnetic field
into the regime where they cannot be neglected.

Here, we run two stages of simulations. For consistency, we start with a WD
structure which allows us to reproduce observations, and consistent equations
of state throughout all phases. The first stage extends the study of Paper I
to a three dimensional star to include large-scale wave modes. The ability of
the magnetic field to suppress instability is a function of the wavelength of
the perturbation. In Paper I, we used a small rectangular domain, which
restricted perturbations to those that can fit across the box, 15km. While we
were able to completely suppress RT, larger wavelengths may still be unstable.
This stage uses 3D magnetohydrodynamics (MHD), with the size of $B$ as a free
parameter, and a simplified burning model, starting from a stage at which
Rayleigh-Taylor instabilities are known to develop. The second stage extends
Paper II by employing the magnetic field topology from a more realistic model,
namely the simulations of the first stage. In the second stage we study the
impact of magnetic field strength and topology on positron escape in the late
stages (1000 days) of light curves and line profiles to constrain the size of
the initial $B$-field. This stage uses 1D radiation hydrodynamics, a more
complete nuclear network of 218 isotopes, time-dependent non-LTE models for
atomic level populations, and 3D photon and positron transport. We assume that
the magnetic topology is set by the deflagration phase, and that the post-
detonation dynamics are entirely determined by the explosion. We then simulate
the light curve as described in Hoeflich et al. (2017b), with the addition of
a three dimensional magnetic field that is passively advected along with the
(homologous) expansion of the supernova. This allows us to more accurately
capture the behavior of the positrons, while remaining computationally
feasible.

Our paper is organized as follows. In Section 2 we will describe the numerical
methods employed, and the setup of the simulations. In section 3, we will
discuss stage 1, results of the early deflagration and the morphology of the
field. In section 4, we present stage 2, and will address the impact of
magnetic fields on late-time spectra and light curves. In section 5, we will
discuss the results and implications for the scenario-dependent observables.
We summarize our results, discuss the possible sources of amplification, and
briefly address the differences between $M_{\rm Ch}$ and helium-triggered
explosions.

Figure 1:  Initial Configuration.  _(left)_ A projection of the initially burned region and magnetic field lines. The black line shows the edge of the star at 1975 km. The initially burned material is distributed in fingers radiating from the center. A total of 80 fingers starting 200 km from the center with length 400 km are used. This initial configuration is identical for all 4 simulations, with only the strength of the field changing. _(right)_ The initial thermal pressure (black line) and magnetic pressures vs. radius for all four simulations. The black line shows the pressure profile (identical for all runs) while the colored power law lines show the magnetic pressure (left axis) or field strength (right axis).  Table 1: Model (run) names vs. the strength of the initial magnetic field, $B$.  Model Name | Magnitude of the Initial
---|---
| Magnetic Dipole Moment
| [$\rm{G\ cm^{-2}}$]
D22 | $10^{22}\hat{x}$
D26 | $10^{26}\hat{x}$
D27 | $10^{27}\hat{x}$
D28 | $10^{28}\hat{x}$

##  2 Method

Our two stages of simulations both begin from the same hydrostatic white
dwarf. The initial model for all simulations in both stages is based on Model
23 from Hoeflich et al. (2017b).

The first stage is restricted to the early deflagration phase; this phase
explores the impact of magnetic fields on the burning, and examines the
evolution of the topology of the magnetic field. This stage follows the non-
distributed regime of burning in three dimensions with 3D magnetohydrodynamics
(MHD).

The second stage also begins from Model 23, and follows the SNIa to late
times, examining the impact of magnetic field strength and topology on the
transport of positrons, and the resulting effect on light curves and spectra.
Stage 2 uses a hybrid approach, employing distributed burning in 1D, while the
3D magnetic field is a passively carried with the explosion. To initiate the
magnetic field in the second stage we map the magnetic field from stage 1 onto
a spherical comoving grid. We follow the evolution to the phase of free
expansion. This two stage approach allows us to use full MHD when the magnetic
field is thought to be dynamically important and to create self-consistent
magnetic topologies, but also follow the light curve for much longer than
possible for an explicit 3D MHD simulation.

The first stage, discussed in Section 2.1.2, uses the 3D MHD code Enzo (Bryan
et al., 2014). For the second stage, discussed in Section 2.2.2, we use our
non-local thermodynamic equilibrium (non-LTE) code for 3D HYDrodynamcal
RAdiation transport HYDRA. Results are discussed in Sections 3 and 4,
respectively.

###  2.1 Stage1: Initial Deflagration

####  2.1.1 White dwarf setup

Each of our four stage 1 simulations begins with a hydrostatic WD with a
central density of $1.15\times 10^{9}\rm{g/cm^{3}}$ and an initial radius of
1975 km. To seed the burning instability, we begin with the central region
initially burned in a star-like pattern. We superimpose a magnetic dipole
oriented along $\widehat{\mathbf{x}}$ at the origin. The only parameter that
changes between the four models is the dipole moment magnitude, i.e. the
strength of the magnetic field. The cartesian domain is
$4000\,\,\mathrm{km}\mathbf{{}^{3}}$, with outflow boundaries on all
quantities.

Model parameters are chosen to resemble the physical conditions of a WD after
the onset of the deflagration stage, while the burning is still in the non-
distributed regime. We interpolate from Model 23 to the higher resolution 3D
grid using Taylor series for the density and the pressure, subject to the
requirement to fit the original profiles and derivatives at the innermost
points.

The internal energy is initialized to maintain hydrostatic equilibrium in the
unmagnetized case. The value of the adiabatic index, $\gamma=1.35$ was chosen
to be consistent with stage 2. This choice is baesd on the equation of state
(EOS) as implemented in HYDRA (see appendix B, and Paper I). The result is a
radial pressure profile, not including magnetic pressure, from which we
construct the internal energy profile. Owing to the combination of the
mismatch between the spherical structure and Cartesian grid and the magnetic
pressure, we were left with small residual velocities within the inner 100 km,
which we discuss further at the end of this section.

Most of the star is pure fuel (50:50 carbon/oxygen). We initialize a burning
front at the center (100% ${}^{56}\rm{Ni}$). To build the burning front, we
first triangulate a sphere with radius 200 km with 80 facets, then add 200 km
fingers on each facet. A projection of this star-like pattern can be seen in
Figure 1.

The initial magnetic field is a global magnetic dipole with a moment along the
$x$-axis. To ensure the divergence of the magnetic field is numerically zero,
we initialize the vector potential and take its curl to produce the magnetic
field. The vector potential is

$\displaystyle\mathbf{A}(\mathbf{r})=\frac{\mathbf{m}\times\mathbf{r}}{4\pi
r^{2}},$ (1)

where $\mathbf{m}=M\hat{x}$. The initial dipole moment $M$ are
$10^{22},10^{26},10^{27},\text{and }10^{28}\,\rm{G}{cm^{-2}}$. We refer to
each simulation by the order of magnitude of the dipole moment, that is
$D22,D26,D27,$ and $D28$, respectively, see Table 1. The initial magnetic
topology can be seen in the left panel of Figure 1. The angle-averaged
magnetic profile can be seen in the right panel of that figure, which shows
the magnetic pressure and gas pressure (left axis) and magnetic field strength
(right axis).

There are two extraneous sources of acceleration in our setup: the first from
errors in mapping the spherical star to the Cartesian grid; and the second
from our non-force-free magnetic configuration. These accelerations are
confined to the inner 100 km, which is well within the initially burned
region, and only imparts a small amount of excess kinetic energy. It should be
noted that simply offsetting the magnetic pressure by reducing the gas
pressure is insufficient, as the magnetic tension
$\mathbf{B}\cdot\nabla\mathbf{B}$ must also be initially zero. This must be
done numerically. While force-free magnetic fields should be used in future
experiments, the perturbation here is not large enough to affect our
conclusions.

####  2.1.2 Stage 1 method: Enzo

For the first stage of simulations, we use the code Enzo (Collins et al.,
2010; Bryan et al., 2014) modified for non-distributed burning. The code is
designed for astrophysical applications supporting a number of physics
processes, including MHD. It solves the Eulerian equations of ideal MHD as
well as the non-distributed burning equations, which are:

$\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\mathbf{v})=0$ (2)

$\frac{\partial\rho\mathbf{v}}{\partial
t}+\nabla\cdot\left(\rho\mathbf{v}\mathbf{v}+\mathbf{I}P-\frac{\mathbf{B}\mathbf{B}}{8\pi}\right)=-\rho\mathbf{g}$
(3)

$\frac{\partial E}{\partial
t}+\nabla\cdot\left[(E+P)\mathbf{v}-\frac{\mathbf{B}(\mathbf{B}\cdot\mathbf{v})}{4\pi}\right]=-\rho\mathbf{v}\cdot\mathbf{g}+\dot{Q}$
(4)

$\frac{\partial\mathbf{B}}{\partial
t}+\nabla\times(\mathbf{v}\times\mathbf{B})=0.$ (5)

Here, $\mathbf{v}\mathbf{v}$ and $\mathbf{B}\mathbf{B}$ are the velocity and
the magnetic field outer products, $\rho$, $\mathbf{g}$, and $\dot{Q}$ are the
density, the gravitational acceleration, and the rate of energy production
from the nuclear burning. Further $E=e+{\rho v^{2}}/{2}+{B^{2}}/{8\pi}$ is the
total energy density, $P=p+{B^{2}}/{8\pi}$ is the total pressure.

The following equation of state closes the system:

$e=\frac{p}{\gamma-1}$ (6)

To speed the calculation and reduce numerical instability, the gravitational
acceleration, $\mathbf{g}$, is computed from the spherically averaged density,

$\displaystyle\mathbf{g}(\mathbf{r})=\widehat{\mathbf{r}}\frac{4\pi
G}{r^{2}}\int_{0}^{r}\rho r^{2}dr.$ (7)

To check the validity of this assumption, we tested this solver against the
FFT-based gravity solver in Enzo, using the strongly magnetized (and least
round) case. Deviations between the two were at most one percent at the outer
boundary, and $<<1\%$ in the interior.

In the non-distributed burning approximation, we assume that the width of the
burning front is small compared to a zone, and that the fuel is burned
instantaneously to product and specific energy, Q. Our burning operator is
based on Khokhlov (1995a). In this model, the flame propagates via diffusion
of the molar burned fraction, $f$

$\frac{\partial f}{\partial t}+\mathbf{v}\cdot\nabla f=K\nabla^{2}f+R,$ (8)

where $f=0$ is pure fuel. $K$ and $R$ are the diffusion rate and reaction
rate, respectively. $R=R_{0}=const$ if $f$ is between the threshold for
burning, $f_{0}=0.3$, and unity, when no burning is possible, and zero
otherwise. The constants are chosen such that the front diffusion speed
$D_{f}=\sqrt{KR_{0}/f_{0}}=100\,\,\mathrm{km/s}$. The energy produced by the
nuclear burning, $\dot{Q}$, needed in 4, is calculated as:

$\dot{Q}=\mathcal{Q}_{\rm{burn}}\frac{d\rho_{\rm{prod}}}{dt}$ (9)

where $\mathcal{Q}_{\rm{burn}}$ is the nuclear energy released per gram fuel
and $\rho_{\rm{prod}}$ is the partial mass density of the burned product. The
relation between $f$ and $\rho_{\rm{burn}}$ can be seen in Equation A7 in
Appendix References.

We use the Constrained Transport MHD module in Enzo(Collins et al., 2010;
Bryan et al., 2014). We have shown that this method preserves
$\nabla\cdot\mathbf{B}=0$ to near machine precision. For the primary evolution
of the MHD equations, we use Li et al. (2008), and the HLLD Riemann solver of
Mignone (2007). To compute the electric field, the Constrained Transport
method of Gardiner & Stone (2005) is used.

###  2.2 Stage 2: Late time evolution

We wish to study the impact of magnetic field strength and topology on the
transport of positrons in late-time ($\sim 1000$ day) supernovae, and their
impact on light curves and spectra. We show that there is a measurable
increase in late-time light curves ($\sim 0.5^{m}$) if the positrons are
locally trapped by the magnetic field, rather than the unmagnetized case which
allows for positrons to more freely diffuse out of the core. This can be used
to estimate lower limits on magnetic field strengths of published supernova,
see Section 4.

####  2.2.1 Remapping and hybrid model

The cost of continuing the Enzo simulations until t=1000 days is
computationally prohibitive; one of these simulations took several days on 64
processors to compute one second of time. As there are $8\times 10^{7}$
seconds in 1000 days, continuing the run would surely extend past the end of
the funding period. Moreover, the physics needed to follow the explosion
beyond the phase of non-distributed burning are not currently available in
Enzo. Much success has been had by employing 1D simulations, which can afford
more complete physics packages (Hoeflich et al., 2017b). However, the
diffusion of positrons is a fundamentally three dimensional process, being
dependent on the strength, topology, and distribution of the magnetic field.
Thus, for stage 2, we adopt a hybrid approach, where the hydrodynamics is
solved in 1D (assuming spherical symmetry), while the positron transport is
3D. The 3D magnetic field is taken as a passive tracer and expands in a
frozen-in manner with the expanding supernova.

The core assumption is that the magnetic field topology is determined during
the subsonic deflagration phase, and the dynamics during the supersonic
detonation phase are determined exclusively by the overwhelming explosion
energy. We show in Section 3 that only the strongest magnetic field has any
impact on the burning front; the three more reasonable cases are quite
spherical and have similar topologies, which justifies our neglect of the
back-reaction of the magnetic field on the burning in Stage 2. While it is
possible that the topology may continue to evolve after the detonation, it is
unlikely that the field will become less tangled, and it is unlikely that the
field will be further amplified. As the decay timescale for a magnetic field
in a fully ionized WD with temperature $T=10^{9}$K is $10^{9}$ yr (see, e.g.,
Choudhuri, 1998), it is unlikely that substantial decay will take place in the
few hundred days examined here. Thus, passive advection of the magnetic field
is sufficient for our purposes here.

The Stage 2 simulations begin from the same spherical WD that the Stage 1
simulations begin with, and makes use of our non-LTE code for HYDrodynamcal
RAdiation transport, HYDRA. Since the mechanism of DDT is not established, it
is treated as a free parameter. The DDT is initiated by mixing of 0.01
$M_{\odot}$ at the burning front. In this simulation, the DDT is triggered
after burning of about $0.27M_{\odot}$ corresponding to a transition density
of $2.5\times 10^{7}\rm{g/cm^{3}}$. Unlike some previous studies (Hoeflich,
2006; Fesen et al., 2007), we will not consider off-center delayed-detonation
transition (Hoeflich, 1990, 1995a, 2003a, 2003b; Diamond et al., 2015; Telesco
et al., 2015; Hoeflich et al., 2017b). This technique produces light curves
and spectra that are consistent with normal-bright SNe Ia (Hoeflich et al.,
2017b). The addition in this current work is the improvement of the positron
transport.

We should note that an attempt was made to directly map the results of Stage 1
to the 3d hydrodynamics models in HYDRA. The structure is close to
hydrostatic, and $P_{gas}$ is many order of magnitude larger than $P_{mag}$
(see Fig. 1). However, omitting even a small magnetic pressure caused
numerical instability on the grid-scale, crashing the simulations. Thus, we
use this more simplified approach.

####  2.2.2 Stage 2 method: HYDRA

For the Stage 2 simulations we use our non-LTE code for HYDrodynamcal
RAdiation transport, HYDRA. Here we briefly outline the modules included. For
more details, see Appendix B.

The simulations begin with a spherical hydrostatic C/O WD, identical to that
of Stage 1 (Hoeflich et al., 2017b). The simulations utilize a nuclear network
of 218 isotopes during the early phases of the explosion; detailed, time-
dependent non-LTE models for atomic level populations; and $\gamma-$ray and
positron transport and radiation-hydrodynamics to calculate low-energy LCs and
spectra (Hoeflich, 1995a, 2003a; Penney & Hoeflich, 2014).

The results use spherical hydrodynamics until 10 days after the explosion.
During the ongoing deflagration phase, the rate of burning is parameterized
based on physical flame models and calibrated to 3D hydrodynamical models by
A. Khokhlov (Domínguez & Hoeflich, 2000; Khokhlov et al., 1997; Gamezo et al.,
2005). The decay of ${}^{56}\rm{Ni}$ increases the specific energy
corresponding to a velocity of $\approx
3{,}000\,\,\mathrm{km}\,\,\mathrm{s^{-1}}$ over the course of its half life,
6.1 days. This ongoing energy input requires the hydrodynamics and the time-
dependent radiation-transport equations to be solved simultaneously in order
to properly treat cooling by the $PdV$ term. For this early phase we solved
the radiation transport equations in co-moving frames using $\approx 100$
frequency groups. The simulation uses atomic models with $\approx 10^{6}$ line
transitions taken from the compilations of Kurucz (1993), supplemented by
forbidden line transitions of iron group elements in the lower ionization
states (Telesco et al., 2015).

After 10 days, the hydrodynamical evolution is assumed to be homologous
expansion, $v\propto r$. This adequate because the kinetic energy dominates
all other energy sources, and the sound speed is many orders of magnitude
smaller than the expansion velocity.

For computational efficiency, we did not recalculate the photosphere and
transitional phase of expansion but skipped to the nebular phase starting at
day 100.

Magnetic fields impose the necessity of 3D treatment for the positron-
transport and, with it, for the 3D radiation transport for both high and low-
energy photons. For Stage 2 simulations discussed here, we use the 3D Monte-
Carlo modules for the photon and positrons.

Because of the assumptions above, we use angular averaged departure
coefficients, i.e. the relative population number of the atomic levels
relative to its ground state within each ion. This is justified because the
infrared (IR) radiation originates from optically thin layers, and
luminosities become isotropic. We assume stationary radiation transport, i.e.
we neglect the implicit time-dependence in the transport terms for photons and
positrons, and in the rate equations. For details, see Appendix B.

We restrict our discussion to day 1000 because after that, additional
radioactive isotopes (e.g. ${}^{57}{Co},^{55}Fe,^{44}{Ti}$) are expected to
dominate the energy input, giving a natural limit for our study.

##  3 Stage 1 evolution: morphology and burning rate during the early
deflagration phase

Figure 2: Slices of Nickle Density for all four simulations. Magnetic field
strength increases downward, and time increases to the right. Only the most
strongly magnetized case shows appreciable modification to the flow
morphology.

Figure 3: Magnetic Field Strength and Streamlines. In our weak-$B$ model, D22
(top row), the magnetic field develops eddies following the Rayleigh-Taylor
(RT) instabilities. The structure of the initial global magnetic dipole is
almost lost. It requires a much higher initial $B$, as in D28 (bottom row),
for it to survive the RT instabilities. Since the latter magnitudes are
unrealistic, this agrees with the lack of observational evidence for
directional dependence, suggested by models in non-DDT scenarios. Figure 4:
Energy vs. Radius for the first $0.5\rm{s}$. Here we show magnetic energy
density (solid lines), gas pressure (black solid lines) and kinetic energy
density (dotted lines) averaged over the sphere. Time increases to the right.
The simulation $D22,\ D26,\ D27$ and $D28$ are color coded as red, yellow,
cyan, and blue, respectively, as in Figure 1. As the pressure profiles are
nearly identical, we plot only one for clarity

###  3.1 Burning and Rayleigh-Taylor

The evolution of the nickel density for each of our four simulations can be
seen in Figure 2. The top row of figures has the weakest magnetic field, while
the bottom has the strongest. Time increases to the right, from the initial
condition to $t=0.6\rm{s}$. In each simulation, the burning proceeds outwards,
retaining an imprint of the initial perturbation. It is notable that the
magnetic field has little impact on the nickel distribution for any but the
most strongly magnetized simulations.

Gas heated by the burning subsequently develops rising fingers of magnetized
gas due to the Rayleigh-Taylor instability. This is seen in the magnetic field
distributions in Figure 3, which shows slices in the most weakly and most
strongly magnetized simulations, $D22$ and $D28$ respectively, at
$t=0.6\rm{s}$. Due to flux freezing, the magnetic energy that is initially
strongest in the center is dragged upwards with the heated “low” density gas.
The instability, and the correlation between the magnetic field and the low
density gas, is most evident in the strongly magnetized run, see the bottom
left panel of Figure 3. This is counter to our initial expectation, which was
that the strongest magnetic field would suppress the Rayleigh-Taylor
instability (Chandrasekhar, 1961; Hristov et al., 2018). Here Rayleigh-Taylor
seems to be more pronounced in the most strongly magnetized run.

This extra RT can be explained in part by the extra energy released during the
settling of the initial conditions. Figure 4 shows the energy densities for
each of the four simulations (with red, yellow, cyan, and blue in order of
decreasing magnetic field strength). Black solid lines show the gas pressure,
color solid lines show magnetic energy, color dotted lines show kinetic
energy. The pressure profile for each can be seen as the black dashed line.
The pressure profile evolves very little during this evolution, and is nearly
the same for each of the simulations. At $t=0\rm{s}$, the magnetic energy is
below the gas pressure for all radii for the two weakly magnetized fields. The
two more strongly magnetized runs have a small excess of magnetic energy. The
excess magnetic energy is releases in the first 0.01s, and in turn drives a
somewhat larger flow than that seen in the other two simulations. It can also
be seen that the strongly magnetized run, $D28$, the magnetic and kinetic
energies are roughly balanced over all radii by $t=0.05\rm{s}$. For the other
simulations, kinetic energy due to the burning outweighs magnetic energy at
outer radii. It is interesting to note that even though the magnetic energy is
large in $D28$, there is still clearly defined Rayleigh-Taylor. This shows
that equipartition of magnetic and kinetic energies is not a sufficient
criterion for suppression of hydrodynamical instabilities.

It can be seen in the bottom row of Figure 2 that the burning Nickel in the
most magnetized run is flattened along the pole-ward direction. This is due
mainly to the anisotropic magnetic pressure, which is larger along the poles.
Only in D28 is this pressure large enough to alter the burning.

Figure 5:  The Auto Correlation Function for the magnetic field alignment.
Black curves show the correlation in the magnetic field direction,
$AC_{\theta,B}$. The velocity direction auto-correlation, $AC_{\theta,v}$, is
shown in red. As expected, correlations in the magnetic field are determined
by correlations in the velocity for all but the most strongly magnetized
simulations. In the most strongly magnetized case, the velocity is more
correlated, and the field is less correlated than the velocity. The blue line
shows $AC_{\theta,B}$ for a simulation of driven turbulence, to indicate an
extreme case of disorder. The correlation in the field itself is intermediate
between an ordered dipole and a fully turbulent field. Figure 6: The auto-
correlation length, defined by eq. (12), of the magnetic field as a function
of time, in units of the initial white dwarf radius, $R_{\rm{WD}}$. Initially,
the magnetic field is correlated over half the star, but decreases to a
fraction of that. Improved resolution would likely cause this to decrease
further.

###  3.2 Tangling

The structure of the magnetic field, while initially a dipole, is quickly
tangled by the kinetic motions of the gas. This can be seen clearly in Figures
3 and 10, which shows the magnetic field lines for the two extreme cases. The
magnetic field structures clearly follow the density structures, becoming
substantially more tangled. This is due to the action of flux freezing,
whereby magnetic flux is pinned to the moving fluid.

In order to quantify the degree of tangling in the magnetic and velocity
fields, we use the auto-correlation functions
$AC_{\theta,\hat{\mathbf{b}}}(\Delta x)$ and
$AC_{\theta,\hat{\mathbf{v}}}(\Delta x)$ of the respective direction vectors,
$\hat{\mathbf{b}}=\mathbf{B}/B$ and $\hat{\mathbf{v}}=\mathbf{v}/v$:

$\displaystyle AC_{\theta,B}(\Delta x)=\frac{1}{4\pi(\Delta
x)^{2}V}\int_{S(\mathbf{x},\Delta x)}ds(\Delta\mathbf{x})\int d^{3}\mathbf{x}\
\hat{\mathbf{b}}(\mathbf{x}+\Delta\mathbf{x})\cdot\hat{\mathbf{b}}(\mathbf{x}),$
(10) $\displaystyle AC_{\theta,v}(\Delta x)=\frac{1}{4\pi(\Delta
x)^{2}V}\int_{S(\mathbf{x},\Delta x)}ds(\Delta\mathbf{x})\int d^{3}\mathbf{x}\
\hat{\mathbf{v}}(\mathbf{x}+\Delta\mathbf{x})\cdot\hat{\mathbf{v}}(\mathbf{x}).$
(11)

These functions measure the degree to which the fields $\hat{\mathbf{v}}$ and
$\hat{\mathbf{b}}$ are correlated with themselves at a separation of
$\Delta\mathbf{x}$. The interior integration is an average over the domain
$V$. The outer integration is over the direction of the shift $\Delta x$, so
the functions depend only on the magnitude of the shift and not its direction.
The above function are normalized so that $AC(0)=1$. At zero separation,
anything is perfectly correlated with itself. A field with more spatial
disorder, and hence increased tangling of the fields, will show a lower
correlation with itself at a certain distance. Figure 5 shows the auto-
correlation function for our dipole (labeled $t=0$) which shows a nearly
linear decrease in the correlation with distance. On the same graph, a
“maximally disordered” configuration of driven high Reynolds number MHD
turbulence (blue line) is shown for comparison. The more tangled turbulent
state shows lower auto-correlation, despite being statistically spatially
homogeneous, due to the chaotic nature of turbulence.

Figure 5 also shows the $B$ and $v$ auto-correlations at $t=1.0\,\rm{s}$ for
all runs. All four simulations begin with the same dipole field. As the
simulation proceeds, the burning drives a number of instabilities, among them
Rayleigh-Taylor and Kelvin-Helmholtz, which cause fluctuation in the field
direction due to flux freezing. This can be seen as a decrease in the auto-
correlation. Three of the simulations have identical auto-correlation
functions. These are shows as dashed, dotted, and dot-dashed lines in black,
but are indistinguishable from each other. Only the most strongly magnetized
run shows any deviation. That simulation (solid black and red lines) shows
small differences in the magnetic correlation, with the magnetic field
resisting some small amount of tangling. It additionally shows a substantial
difference in the velocity correlation, with larger scale flows being more
active.

We can further quantify the tangling by defining the auto-correlation length,

$\displaystyle L_{ac}=\int_{0}^{\infty}\mathbf{d}\Delta xAC_{\theta,B}(\Delta
x)/AC_{\theta,B}(0),$ (12)

which describes the length at which there is appreciable correlation. This
quantity can be seen relative to the initial white dwarf radius in Figure 6.
This plot shows that the correlation length monotonically decreases with time,
as the field becomes more tangled, and this length is not particularly
sensitive to our initial magnetic field.

The conclusion from this discussion is that the topology of the magnetic field
is determined primarily by the burning front, and not greatly influenced by
the initial field. The topology of the magnetic field, as well as its
strength, influences the rate at which positrons diffuse away from their place
of origin. Larger correlation lengths allow positrons to travel further, while
a more tangled field traps the positrons. As our initial conditions are
heavily idealized, we cannot predict what values $L_{ac}$ will attain for a
real SNIa, but we have established that a burning supernova cannot support a
magnetic field topology as simple as a dipole.

Additionally, we are further justified in neglecting the back-reaction of the
magnetic field on the flame during the Stage 2 simulations. Since the strength
of the field does not impact its topology, it also does not impact the burning
front in an appreciable way, and we can neglect it for our purposes in Stage
2. While there is a small effect of the field in the strongly magnetized run,
it is unlikely to change our conclusions in a qualitative manner.

Figure 7: Nickel mass vs. time for all four simulations _(top)_ and mass
difference relative to the D22. The three weakly magnetized simulations are
nearly indistinguishable, but the strongly magnetized run shows a sustained
reduction of the burning rate for the duration of the simulation.

###  3.3 Burning

We show the difference in burned mass in Figure 7. From this, we see the
strongly magnetized case, $D28$, burns about 15% less nickel than the other
three runs. For the initial drop before 0.1s, the difference in burning rate
is due to the initial excess of magnetic energy that decreases the central
density somewhat. Beyond 0.1s, the difference in burning rates is largely due
to differences in the morphology of the interface. This can be seen by
comparing the difference in burning rate (the slope of the bottom panel of
Figure 7) and the difference in $L_{\rm{AC}}$ in Figure 6.

Figure 8: Magnetic energy vs. time for each of the runs _(top)_. _(Bottom)_
Relative magnetic energy change after the initial conditions relax. No
evidence of field growth is seen except for a small increase in $D27$.

###  3.4 Amplification

Figure 8 shows the magnetic energy vs. time for each simulation. There is an
initial drop as the field relaxes from our out-of-equilibrium initial
conditions, and the field energy is converted in part to kinetic energy (see
Figure 4). Since the burning is nearly identical, the difference in kinetic
energy between the runs comes from the release of magnetic energy. To examine
the subsequent evolution, we plot the change in total magnetic energy relative
to $t=0.1$s, right after the initial conditions relax. The total magnetic
energy for three of the runs continues to decrease, while the $D27$, the
second-strongest field, shows a slight increase.

Figure 9: Density vs. Radius at $t=0.6\rm{s}$ for each of our simulations.
Error bars show deviation from spherical. _(Bottom)_ The standard deviation in
the density. Errors in the strongly magnetized run are on the order of a few
percent.

###  3.5 Spherical Symmetry

Even though there is substantial magnetic field, the structures within the WD
stay relatively spherically symmetric on average because the overall system
stays in hydrostatic equilibrium dominated by the gravitation. Figure 9 shows
the angle-averaged density vs. radius at $t=0.6$ s. The top panel shows the
average density and the standard deviation from spherical,
$\sigma_{\rho}(r)=\sqrt{\langle(\rho(r)-\bar{\rho}(r))^{2}\rangle}$. Error
bars are shown on all plots, but smaller than the line width for small radius.
The deviation from spherical is most pronounced in the strongly magnetized
run. The bottom panel of that figure shows the normalized standard deviation,
$\sigma_{\rho}/\rho$. The peak $\sigma_{\rho}$ is $8\%$ for the most
magnetized case. The few percent deviation enjoyed by all simulations at
$r>3\times 10^{7}\rm{cm}$ is due to the boundary of the domain. Future work
will employ a slightly larger domain to avoid this effect. The majority of the
simulations are spherical to the fraction of a percent for the bulk of the
evolution.

Figure 10: Magnetic field line “hairballs” In our weak-$B$ model, D22 (left
column), the velocity field develops Rayleigh-Taylor (RT) instabilities and
the magnetic field eddies follow. The structure of the initial global magnetic
dipole is almost lost. It requires a much higher initial $B$, as in D28 (right
column), for it to survive the RT instabilities. Since the latter magnitudes
are unrealistic, this show why there is no observational evidence of
directional dependence, suggested by models in non-DDT scenarios. For clarity
the figures show field lines originating in the $x<0,\ y>0,\ z>0$ octant only.

Figure 11: Light Curves. Influence of the magnetic field morphology on light
curves for large-scale undisturbed dipoles (_dipole_) and the field structure
produced by the turbulent morphology produced by Rayleigh-Taylor instabilities
starting from and initial dipole field in the WD for model D22 (solid line)
and D28 (dotted). The left plots shows the light curve for a fixed dipole that
does not evolve (solid) and the turbulent field taken from D22. The right plot
shows the evolution of the deviation with respect to our fiducial light curve
to show the differences between D22 (solid line), D28 (dotted line) and the
fixed dipole (dashed). As background explosion model, we use the delayed
detonation model 23 with parameters from Hoeflich et al. (2017b) for a normal-
bright Type Ia Supernovae, and imprinted the B-morphology with surface fields
up to $10^{13}G$ which causes local trapping of positrons (see Table 2). We
show the corresponding bolometric light curves (LCs) (left), and the
difference (right) relative to the case of full local trapping. Note that a
field size of ${10}^{11}G$ is already close to local trapping even at 1000
days. The details of the morphology of a turbulent field has minor impact (D22
vs. D28), but large scale ordered fields allow for a significantly larger
fraction of positrons to escape. As a result, the solution for the dipole
field of $10^{3}G$ is almost indistinguishable from the field-free case and
its field of $10^{6}G$ produces an escape similar to a turbulent field but
some three orders of magnitude smaller. In addition, we give the loci of the
upper limits for SN2012cg, SN2003hv, SN2014J, SN1992A and SN2011fe (blue dots
from left to right). All SNe are consistent with complete trapping of
positrons. Note that an assumed uncertainty of $0.15$ to $0.2^{m}$ for
bolometric LC implies the need for observations beyond $\approx 500days$. As
discussed in the text, all LC observations follow the case of local trapping
of positrons (green dotted). At about 700 to 1000 days, additional radioactive
heating places their locus above the full trapping case limiting the use of
LCs to times before day 1000. For the estimates of the lower limit in $B$
(Table 3) we assume data starting beyond $\approx 250$ days using $L_{bol}$
for a normal bright SNe Ia (Paper II). Note beyond day 300 more than 95% of
the energy input is from positron (see text).

##  4 Stage 2 Evolution: Observable consequences and implications of late time
evolution

Based on the results described in the previous section, we will show the MHD
effects discussed above based on the imprint on the late-time light curves
(Milne et al., 2001; Penney & Hoeflich, 2014) and NIR line profiles (Hoeflich
et al., 2004; Penney & Hoeflich, 2014; Diamond et al., 2015) in the framework
of delayed-detonation models (Khokhlov, 1989; Hoeflich & Khokhlov, 1996).

In the stage 1 simulations, we studied the evolution of the morphology of the
$B$-field based on the initial density structure and the adiabatic coefficient
$\gamma$ consistent with the equation of state used in the explosion
simulations of HYDRA, namely the parameters of the spherical delayed-
detonation model 23 of Hoeflich et al. (2017b). In stage 1, the size of the B
is a free parameter. This model starts with a WD originating from a main
sequence star of $7~{}M_{\odot}$ and solar metallicity $Z_{o}$. The accretion
rate from the donor star has been tuned to produce a thermonuclear runaway at
a central density $\rho_{c}=2\times 10^{9}g/cm^{3}$. In stage 2, we followed
the deflagration phase during the distributed regime of burning and the
subsequent detonation phase. The delayed-detonation transition has at a
transition density $\rho_{tr}=2.3~{}10^{7}~{}g/cm^{3}$. We triggered the
delayed-detonation transition ’by hand’ because the mechanism(s) leading to
the DDT have not been established. We choose this model because it can
reproduce the observed curves, color-magnitude relation, optical to mid-IR
spectra of a typical normal-bright SNe Ia, namely SN2014J (Diamond et al.,
2015; Telesco et al., 2015; Hoeflich et al., 2017b).

As in previous studies on the effects of $B$ on positron and photon transport,
we consider the nebular phase when the optical depth is and densities are low,
wherein forbidden atomic transitions dominate. The instant energy input from
radiative decay of ${}^{56}Co$ governs the bolometric luminosity. Moreover,
during the late phase, the envelope is homologeously expanding, so that
position and velocity of a mass element are related as
$r(m,t)\propto|\vec{v}(m,x,y,z)|\times t=v(m)\times t$. Only the radial
component of $\vec{v}$ remains. The line profile of unblended features is a
result of the Doppler shift of the mass element $m$ and the energy deposition.

The HYDRA code employs full-3D transport and hydrodynamics including the time
dependence (see Appendix B). However, the results in the nebular phase are
dominated by forbidden lines in an mostly optically thin envelope. The
solution of the statistical equations is dominated by the branching ratio
between atomic energy levels $j$ and $i$, here determined predominately by the
Einstein values $A_{ji}$ (Motohara et al., 2006; Penney & Hoeflich, 2014;
Diamond et al., 2015; Telesco et al., 2015), and the radiation transport is
close to the optically thin limit leading to the stability of the results.

As discussed below, the _lower_ limits given below for the size of $B$-field
derived from spherical LCs is not affected by ${}^{56}Ni$ plumes or
instabilities. This is because mixing out of ${}^{56}Ni$ would increase the
escape probability of positrons and $\gamma$-rays (see Fig. 11) which would
increase the estimate for $B$, but our lower limits will still be valid. The
lower limits correspond to average surface fields assuming a magnetic topology
identical to $D22$, for which our simulations of stage 1 show small deviations
from sphericity. Even using the slightly less spherical model $D28$ results in
only a slightly higher field estimate, see Figure 11.

Note that we make the case that asphericities in $\rho$ for the highest
$B$-fields may not survive the further evolution based on pure 3D hydro-
simulations (Khokhlov, 2001; Gamezo et al., 2003b, 2005; Röpke et al., 2012).
However, this has not been shown in this work.

We use the morphology of the magnetic field for model D22 at 0.6 seconds, and
scale the average surface magnetic field between 1 and $10^{13}G$. For a grid
of magnetic fields of $10^{0,1,3,4,6,9,11,13}$ G, the positron, $\gamma$, and
photon transport problem has been simulated at 100, 150, 200, 300, 400, 500,
750 and 1000 days. For verification of our main conclusion, we have used the
imprint of the highest magnetic field hydro-model (D28) with $10^{13}G$. We
note that large scale asymmetries in the density and chemical structure will
affect the line profiles, as shown in previous simulations for 3D turbulent or
off-center delayed-detonation models (Höflich, 2002; Höflich et al., 2006;
Motohara et al., 2006).

###  4.1 Magnetic field effects on late time light curves

The results can be well understood within the following framework. $B$-fields
are frozen in comoving frame because of the plasma. Positrons are scattered
and lose energy by collisions with electrons till the energy is below the
binding energy. At this point, they are annihilated. For details, see appendix
B.

To first order, the SN envelope is freely expanding. The density drops with
time $t$ as ${1/t^{3}}$ and the length scale increase $\propto t$. As a
result, the mean free path of positrons relative to the expanding envelope
increases like $t^{2}$ before annihilation.

In the presence of $B$-fields, positrons gyro around the field lines with the
Larmor radius. If the structures in the field morphology are smaller than the
Larmor radius, positrons are trapped. Because structures grow linearly with
$t$ but the Larmor radius increases like $t^{2}$, positrons will eventually
propagate beyond these structure (see Table 2).

Table 2: Maximum Larmor radius as a fraction of the size of the supernova envelope at various times and initial surface field strengths. During the nebular phase, the expansion is homologous with the distance of a mass element $r(m)\propto v(m)\times t$. $v(m)$ and $r(m)$ are equivalent measures of the envelope structure. However, $v(m)$ corresponds to a spectral Doppler width observable in spectral line profiles. Here, the size of the envelope is defined by the layer expanding with a velocity of 40,000 km/sec. In our simulations, the typical size of the $B$-field structure is produced by RT instabilities ($\approx 500...1,000km/sec$ in velocity units or $1...2\times 10^{-2}$ in units of the envelope size.) field/time | 100d | 300d | 500 d | 750 d | 1000d
---|---|---|---|---|---
$10^{3}G$ | 1.9 | 5.7 | 9.6 | 14.4 | 19.1
$10^{4}$G | .19 | .57 | .96 | 1.4 | 1.9
$10^{6}$G | .001 | .005 | .009 | 0.013 | 0.02
$10^{9}$G | 1.9E-6 | 5.E-6 | 9.E-6 | 1.3E-5 | 2.E-5
$10^{13}$G | 5.E-11 | 1.2E-10 | 2.E-10 | 3.3E-10 | 5.E-10

Theoretical bolometric LCs, $L_{bol}$, and the difference in magnitudes in
$V$, a proxy for $L_{bol}$, are shown in Fig. 11 for magnetic fields between 0
and $10^{9}...^{11}$G. The zero-field case is identical to the
$B=10^{3...4}\rm{G}$ dipole field, and the total local trapping case is
virtually identical to the $B\geq 10^{9}\rm{G}$ turbulent cases for all times.
The early LCs are identical because $\gamma-$ray transport does not depend on
$B$ and the free mean path of positrons is small. At about 220 days, the
energy input from positrons equals that from $\gamma$-rays. Light curves with
different magnetic fields start to diverge after about day 300 when the mean
free path some of the positrons becomes large compared to the envelope. Small
scale structures and high $B$-fields trap positrons for longer than low $B$
and or large scale dipoles. The luminosity for the simulation with $B=10^{9}G$
and a turbulent field follows closely the energy input from the decay of
${}^{56}\rm{Co}$ even at day 1000, which shows that the positrons are not able
to escape the region of their birth. The difference (Figure 11, right) between
the models may be in excess of $0.5^{m}$ which should be easily measurable
given current observational accuracy. We examine observed LCs in the next
section.

####  4.1.1 Magnetic field strengths in known supernovae

We want to combine the results of Figure 11 with observations. Using our
simulations (Fig. 11, right plot) and applying the deviation of the observed
range in time, a given SNe Ia provides estimates for the size of the
$B$-field. The detailed reconstruction of $L_{bol}$ is beyond the scope of
this paper and, thus, we refer to $L_{bol}$ (or $V$ as proxy) given in
literature.

We use the radioactive tail of the fiducial light curve that keeps the
positrons local, and compare it to simulations with a variety of magnetic
field strengths and topologies. We then compare these differentials to
observed supernovae to determine lower limits of their field strengths. This
can be seen in the right panel of Figure 11, which shows the deviations from
full trapping for each of our model simulations, as well as the upper limits
on the deviation (thus lower limits on the field strength) for several known
supernovae.

Relatively few late-time LCs are available. Milne et al. (2001) suggested some
indication for positron escape by day $\approx 170$, but this was based on
strictly radial $B$-fields in the models, and neglected near and mid-IR
emission. Multi-band LCs are required to reconstruct $L_{bol}$ with sufficient
accuracy to include the redistribution of photons to the NIR and mid infrared
(MIR) within an accuracy of $\approx 0.15^{m}$ (Stritzinger et al., 2002;
Sollerman et al., 2004; Gall et al., 2018). If NIR data are not available, we
use the redistribution functions from our models for normal bright and
transitional SNe (Hoeflich et al., 2017b; Gall et al., 2018) and
redistribution functions based on observations of SN 2014J (Telesco et al.,
2015). We use the intrinsic uncertainty from the scatter in the measurements.
For $L_{bol}$, based on the references above, we add to the measurement error
in the LC points an additional error of $0.15^{m}$ to the reconstruction of
$L_{bol}$ if IR observations are available, and $0.2^{m}$ if they are not.

Table 3: Lower limits for magnetic field strengths from known supernovae. The dates correspond to the last data point observed. The loci in Fig. 11 are the last data point before additional energy sources contribute to the light curve.  Name | t [days] | $B_{\rm{turb}}$ | $B_{\rm{dipole}}$ | reference
---|---|---|---|---
SN1992A | $926$ | $10^{5.5}$ | $10^{8.5}$ | Cappellaro et al. (1997)
SN2003hv | $800$ | $10^{6}$ | $10^{9}$ | Leloudas et al. (2009)
SN2011fe | $930$ | $10^{5}$ | $10^{8}$ | Kerzendorf et al. (2014)
SN2012cg | $640$ | $10^{4}$ | $10^{7}$ | Graur et al. (2016)
SN2014J | $800$ | $10^{5.5}$ | $10^{8.5}$ | Yang et al. (2018)

None of the observations fall below the light curve expected from the decay of
${}^{56}\rm{Co}$, indicating full local trapping by magnetic structures. To
determine the lower limits of the initial $B$, we use the curves of Figure 11
and determine the maximum deviation from the theoretical curves which are
consistent with the decay of ${}^{56}\rm{Co}$. Obviously, this limits depend
on the coverage of the late-time LCs and the error bars. Observed LCs include
SN1992A ($\approx 950d$, Cappellaro et al. (1997)), SN2003hv ($\approx 700d$,
Graur et al. (2016)), SN2011fe ($\approx 930d$, Kerzendorf et al. (2014)), and
SN2012cg ($\approx 1060d$, Graur et al. (2016)). SN2014J has been observed
even beyond 1200 days, showing the onset of radioactive decay of elements
besides ${}^{56}\rm{Co}$. This flattens the LC by about day 800, indicating
the onset of ${}^{57}Co$ decay as the dominant energy source (Yang et al.,
2018). SN2014J was found to be consistent ${}^{56}\rm{Co}$ and with no
signature for positron escape nor an IR catastrophe predicted by Fransson &
Sonneborn (1994). Similarly, SN2003hv, SN2011fe are observed after onset of
${}^{57}\rm{Co}$ heating or interaction, and neither LC drops below the
${}^{56}\rm{Co}$ or ${}^{57}\rm{Co}$ line, giving evidence for full positron
trapping. Results for the lower limit are summarized in Table 3. This shows
that high $B$-fields are common, well in excess of typical fields observed in
WDs (see Sect. 1). Note that, in the presence of a deflagration phase in
$M_{\rm Ch}$-mass explosions, large scale dipole fields will be transformed
into small-scales turbulent fields. Alternative explosion scenarios involve
pure detonation burning. For these, the dipole field morphology may remain
(see next section) leading to even larger lower limits for the initial $B$
field.

Figure 12:  Angle averaged line profiles for the ’almost’ unblended feature
produced by the forbidden [Fe II] at $1.644\mu m$ for an initial average
surface field of $10^{6}G$. The time evolution for the dipole (top row) and
turbulent field (bottom row) shows little evolution before 400 days (left
column) due to the local trapping of positrons. For dipole and turbulent
fields, the profiles start to change after about 400 and 700 days,
respectively (right column). A total of $\approx$ 100,000 line profiles have
been calculated (magnetic field*morphology*phases*$\theta$ and $\phi$
direction $\approx 10^{5}$), are available on request, and will be part of a
data base for different explosion models.

Figure 13:  Same as Figure 12 but [Fe II] line profiles for the dipole field
as seen from $0^{o}$ (_1_), $30^{o}$ (_2_) or from the pole (_3_) at given
times (left). The late-time differences in the profiles ($\sim
1,000...15,000~{}km/s$), although not prominent on this plot, are significant,
given a typical spectroscopy resolution of $\approx 100~{}km/s$. In addition
and as proof of concept, differences between NIR spectra between the dipole
seen equator on and turbulent field are shown on the right with spectra
normalized to the average spectra at day 500. Dipole fields produce a
characteristic, large scale pattern and signals of about $\pm 5\%$ between
$\pm 15,000km/sec$. For turbulent fields, variations correspond to the mass of
the individual plumes and show up as wiggles on the 1% level. Note that
spectroscopy allows precise measurements of Doppler shifts but flux
differences on the 1% level will be reserved to observations with the upcoming
JWST. For details, see text.

####  4.1.2 Magnetic field effects on line profiles

LCs are a powerful tool but, their application is limited as tools to
constrain $B$ because LCs tend to flatten after 700 to 1000 days. Line
profiles and their variations as a function of time have been suggested as an
alternative diagnostic tool in Paper II.

In particular, [FeII] at $1.644\mu m$ has been identified having only minor
blends in the line wing (Hoeflich et al., 2004; Motohara et al., 2006; Penney
& Hoeflich, 2014; Diamond et al., 2015, 2018). All other strong features at
shorter wavelengths are blends produced by multiple transitions. Alternatively
to the NIR [FeII], the MIR has been demonstrated as source for unblended [Co
II] and [FeII] lines by (Telesco et al., 2015) but, likely, those are beyond
reach but for the upcoming JWST. Here, we want to discuss the NIR [Fe II]
(Figures 12 and 13).

The profile is a measure of the distribution of the energy input by positrons
convolved with the Fe abundance. The line profiles are more sensitive to
positron transport effects because changes do not require the escape of
positrons, but non-locality of the positron decay is sufficient (Penney &
Hoeflich, 2014). It provides information about the distribution and possible
off-set of radioactive ${}^{56}Ni\rightarrow^{56}Co\rightarrow^{56}Fe$ and, at
ultra-late times, other decay channels of radioactive isotopes of the iron
group (Hoeflich et al., 2017c). For the time series, we use the small
variation in the line profile as indicator for $B$-fields.

Inherently, the effects of $B$ can be detected earlier, namely by day 400…500
and, in principle, without restrictions to late times well beyond 1000 days.
However, optical depth effects, namely the underlying optically thick
photosphere and non-local energy input by $\gamma-$ rays for iron close to the
center, leads to a change of the forbidden [Fe II] line profile at 1.644 $\mu
m$, rather than revealing positron transport effects (Fig. 12, left). However,
the line profile hardly changes between 300 and 400 days. This requires a
reference point for the spectra later than 300 days. Note that the two
features to the right and left correspond to [Co III] and [FeII] blends and
variations that are caused by the nuclear decay with time (Höflich et al.,
2004). In practice, one problem is related to the costs of obtaining NIR
spectra of high quality, and the need for a time series, which limits the
number of objects to a few (see Sect. 1).

As an example and for a field of $B=10^{6}G$, the evolution of the profiles is
shown as a function of time and orientation in Fig. 12. For both turbulent and
dipole fields, the [FeII] profiles remain almost unchanged between 300 and 400
days after the explosion, because $\gamma-$rays do not contribute
substantially to the energy input, and positrons annihilate locally. After
about day 500, the profile becomes increasingly narrow because positrons at
the outer, high-velocity, low-density regions of the ${}^{56}Ni$ distribution
become non-local, that is they move away from their generation region and
excite other lines and elements within the envelope.

One fundamental problem to decipher the distribution of the radioactive
${}^{56}Ni$ is that a large scale $B$ field, e.g. a dipole, causes strong
directional dependence when a particular SNe Ia is observed. The dependence of
the line profile on viewing angle can be seen in Figure 13. In our example of
a large-scale dipole, positrons can travel almost unhindered along the axis of
symmetry, the polar direction, thus the effective mean-free path has a small
dependence on $B$ as opposed to the highly inhibited travel along the equator.
The changes between 500 and 750d and 500 to 1000d are comparable when seen
from the equator and pole, respectively. The rate of change does depend on
both the size of $B$ and the direction of the observer.

For the turbulent field, the profiles show little directional dependence.
Fluctuations are restricted to the scale of the convective eddies which, in
principle, may be recovered by high signal-to-noise observations with a
resolution of the characteristic eddy-velocities, i.e. 1000 km/sec. We note
that the evolution of the profile, e.g. the mean-half width is non-monotonic
(Diamond et al., 2015) and, thus, requires at least 3 observations to recover.

One of the main result of our MHD simulations is that for a wide range of $B$
the morphology is given by a turbulent field even for a large-scale initial
dipole during the deflagration phase. This does not rule out large dipole
fields in, e.g., pure detonation scenarios.

##  5 Final Discussion and Conclusions

We have studied the effects of magnetic fields on deflagration fronts and
late-time spectra in Type Ia supernovae. We find that small scale magnetic
fields develop as a result of the motions of the burning front, and that
observations put the magnitude of magnetic fields into a regime which may be
both relevant and larger than commonly observed in single WDs.

We employed initial maximum fields between $B=10^{7}\rm{G}$ and
$10^{15}\rm{G}$, in the framework of delayed detonation models. We examine
their imprint on light curves and spectral properties during the nebular
phase. Using our model light curves, we estimate lower limits of
$10^{6}\rm{G}$ for several published supernovae. We also find that Rayleigh
Taylor instabilities develop early, even in the presence of magnetic fields
that, energetically, should suppress the instability.

The two weakly magnetized runs, D22 and D26, are almost entirely identical in
their evolution in all measurable quantities. This is unsurprising, as the
magnetic pressure is several orders of magnitude lower than the gas pressure
at all radii for both of these simulations. This also shows that we are not
troubled by systematic variations other than those caused by the magnetic
field. Tangling develops early (0.5 s) in all runs. There is slightly less
tangling in the strongly magnetized simulations, as well as more pronounced RT
instability. The strong initial field in D28 supports large-scale flow
structures, which results in faster increase in rising plumes as well as
reduced burning surface area. This can explain why D28 has 15% less
${}^{56}{\rm{Ni}}$ vs. lower $B$ as in Hristov et al. (2018).

We do not see any evidence of magnetic field amplification in these
simulations as expected from previous studies (Hristov et al., 2018), though
this does not rule out amplification with other initial magnetic topologies.

We started our simulations with a burned central region which, in size,
corresponds to about 1 to 2 seconds after the thermonuclear runaway in models
of central ignition. This is the onset of RT in models with central ignition
(Khokhlov, 2000; Gamezo et al., 2005). In off-center and multi-spot ignitions,
the development of RT instabilities is faster (Hillebrandt et al., 2013).
Though we do not follow the deflagration front during the regime of
distributed burning, the tangling cannot be expected to decrease. Rather, in
the absence of a large scale flow to order the field, it can only be expected
to increase during the regime of distributed and detonation burning. We show
that the density distribution remains spherical to within $<1\%$ in the most
strongly magnetized case. This is due to the dominance of pressure equilibrium
and gravity, a result found previously by hydro-simulations for quasi nuclear
statistical equilibrium (QSE) of explosive oxygen burning regions (Gamezo et
al., 2003b, 2005).

In Section 4, light curves between 500 and 1000 days are shown to be able to
constrain the strength of the magnetic field, and we use these model light
curves to estimate magnetic field strengths of $>10^{6}\rm{G}$ from several
observed supernovae. Line profiles, or rather their lack of time-evolution,
can be used to constrain $B$ starting at about day 300. For magnetic fields
larger than $10^{6}\rm{G}$, positron are trapped within the magnetic field up
to about day 500 which presents the current limit from observations. However,
changes in the profiles or the width of the line can be used to determine the
magnetic field strength for larger initial fields because, in principle, the
method of using specific line profiles can be applied to even later phases
than light curves. Currently, only a few SNe Ia with sufficiently late time
NIR spectra or with optical and IR light curves are available. A larger number
is needed to confirm whether high $B$ fields are generally present. Figure 11
may provide a tool to get a $1^{st}$ order estimate $B$ fields of future
observations.

As mentioned in the introduction, an obvious question is the origin of such a
field. The magnetic decay timescale of $10^{9}\rm{K}$ plasma is
$10^{9}\rm{yr}$, so it is natural for such a plasma to support magnetic
fields, and the question is identifying its source. Within dynamo theory (e.g.
Brandenburg & Subramanian, 2005) there are several mechanisms to amplify the
field up to the saturation strength of $\sim 10^{14}\rm{G}$ (Chandrasekhar,
1956; Chandrasekhar & Prendergast, 1956; Mestel, 1956; Hristov et al., 2018),
and the nature of the amplification mechanism will be imprinted in the
magnetic field structures. This may happen at several points during the
evolution of the WD. In the Chandrasekhar mass, $M_{\rm Ch}$, explosions
considered here, a WD close to equilibrium begins to burn as a result of
compressional heating, which in turn results from accretion from a companion.
This leads to subsonic deflagration, which then transitions to supersonic
detonation. The correlation length of the flow during each phase will set the
correlation length of the magnetic field, which in turn impacts the escape of
positrons. During the accretion phase, the dominant length scale is the radius
of the white dwarf. During the late-stage run up to the deflagration, the so-
called _smoldering phase_ (Hoeflich & Stein, 2002; Zingale et al., 2011) the
dominant length scale is the pressure scale-height of the star. During the
supersonic explosion, the scale is set by the sound-crossing time scale,
because the flame propagates as a weak detonation. Our results show that
amplification during the deflagration phase is unlikely. Hydrodynamical
simulations (Gamezo et al., 2003b; Röpke et al., 2012) have shown that the
instabilities that could give rise to a dynamo are frozen out during the
expansion phase of the WD. This leaves the accretion phase or smoldering phase
as the likely candidates for magnetic field amplification in WDs.

This leads us to the discussion of these results in the framework of
alternative explosion scenarios. Likely, a combination of different explosion
scenarios are realized in nature. The quest for the predominant mechanism is
still under debate with time changing favorites. This is not too surprising ,
because it is nuclear physics that determines the progenitor structure, the
explosion and the subsequent LCs and spectra (Hoeflich et al., 2013) and
combined with an almost overlapping mass range for the progenitor WD among
different scenarios. This phenomenon is also known as ‘stellar amnesia’. The
mass ranges for normal bright SNe Ia, helium detonations are $\approx
1.0~{}...~{}1.1M_{\odot}$ (Pakmor et al., 2012b; Shen et al., 2018), and 1.28
to 1.38 $M_{\odot}$ for $M_{\rm Ch}$-mass explosions (Hoeflich et al., 1998;
Diamond et al., 2015). For a longer review of explosion scenarios, details and
further references are given in the introduction, Hoeflich (2017) and, for
individual aspects of the explosion, in the subsequent sections under the
chapter titled ”Physics of Thermonuclear Supernovae” in Alsabti & Murdin
(2017). During dynamical merging of two WDs (e.g. (Pakmor et al., 2011b)),
high $B$ fields are likely to develop depending on details of the dynamical
merging process. However, no significant $B$ field amplification can be
expected for He-tiggered explosions because they explode promptly on time-
scales of a second triggered by a supersonic shock, and the progenitor
evolution involves two WD without a deflagration or smoldering phase. However,
a large scale dipole field may develop during the accretion phase similar to
massive He/Co systems described in Pakmor et al. (2021).

Finally, we would like to discuss some limitations beyond those just mentioned
which will be overcome in future. Our study is not a full end-to-end
simulation of the explosion of a WD. Rather, we studied the change of the
morphology of an initial magnetic field during the regime of non-distributed
burning, and the observable consequences. The initial condition of the 3D
simulation starts from a WD structure but a rather advanced stage of burning
to shorten the time till the RT-instabilities can develop. We do not follow
the ignition process. At the end of the 3D simulations, the further evolution
is followed by spherical radiation-hydro without magneto-hydrodynamical
effects.

For our analysis of the NIR-spectra and light curves we have considered a
specific explosion model and realization. A more comprehensive study including
spectra at 750d and beyond is under way and encompassing transitional and sub-
luminous SNe Ia and alternative explosion scenarios.

NIR spectra of SN2004du, 2005df and 2014J, and optical LCs of 5 SNeIa
discussed in this paper suggest $B$ fields in excess of $10^{4...6}G$ to be
common. One obvious limitation is the lack of observations for a larger
sample. With the upcoming James Webb Space Telescope, unblended line profiles
of Co in the MIR will become available (e.g. Gerardy et al. (2007); Telesco et
al. (2015)), and WFIRST will cover the LC evolution of many local SNe Ia ’by
chance’.

We thank the referee for carefully reading the manuscript and many helpful
suggestions. We acknowledge the support by National Science Foundation (NSF)
grant AST-1715133 including the salary for a postdoc. Most of the development
and simulations have been done on the local cluster of the FSU astro-group
including the data storage. This work used the Extreme Science and Engineering
Discovery Environment (XSEDE, https://www.tacc.utexas.edu, Towns et al.
(2014)), which is supported by the NSF grant ACI-1548562, under the XSEDE
allocation TG-AST140008.
_Software_ : Enzo (Bryan et al., 2014), HYDRA (Hoeflich, 1990; Hoeflich et
al., 1993; Höflich, 2003a, b, 2009; Penney & Hoeflich, 2014; Hoeflich et al.,
2017a). Plotting and analysis was done using yt (Turk et al., 2011),
matplotlib (Hunter, 2007), and numpy (Van Der Walt et al., 2011).

## References

  * Adams et al. (1971) Adams, T. F., Hummer, D. G., & Rybicki, G. B. 1971, J. Quant. Spec. Radiat. Transf., 11, 1365
  * Alsabti & Murdin (2017) Alsabti, A. W., & Murdin, P. 2017, Handbook of Supernovae
  * Ambwani & Sutherland (1988) Ambwani, K., & Sutherland, P. 1988, ApJ, 325, 820
  * Athay (1972) Athay, R. G. 1972, Radiation Transport in Spectral Lines
  * Avrett & Loeser (1988) Avrett, E. H., & Loeser, R. 1988, ApJ, 331, 211
  * Benetti et al. (2005) Benetti, S., Cappellaro, E., Mazzali, P. A., et al. 2005, ApJ, 623, 1011
  * Benz et al. (1990) Benz, W., Cameron, A. G. W., Press, W. H., & Bowers, R. L. 1990, ApJ, 348, 647
  * Berger et al. (1998) Berger, M. B., Hubbell, J. H., Seltzer, S., Coursey, J., & Zucker, D. 1998, NBSIR, 87-3597, ”http://physics.nist.gov/PhysRefData/Xcom/Text/XCOM.html”
  * Branch et al. (2005) Branch, D., Baron, E., Hall, N., Melakayil, M., & Parrent, J. 2005, PASP, 117, 545
  * Branch et al. (2009) Branch, D., Chau Dang, L., & Baron, E. 2009, PASP, 121, 238
  * Brandenburg & Subramanian (2005) Brandenburg, A., & Subramanian, K. 2005, Phys. Rep., 417, 1
  * Brooker et al. (2021) Brooker, E., Plewa, T., & Fenn, D. 2021, MNRAS, 501, L23
  * Bryan et al. (2014) Bryan, G. L., Norman, M. L., O’Shea, B. W., et al. 2014, ApJS, 211, 19
  * Buchler (1979) Buchler, J. R. 1979, JQSRT, 22, 293
  * Buchler (1983) Buchler, J. R. 1983, J. Quant. Spec. Radiat. Transf., 30, 395
  * Cannon (1973) Cannon, C. J. 1973, JQSRT, 13, 627
  * Cappellaro et al. (1997) Cappellaro, E., Mazzali, P. A., Benetti, S., et al. 1997, A&A, 328, 203
  * Castor (2004) Castor, J. 2004, Radiation Hydrodynamics (Cambridge: Cambridge Univ. Press)
  * Castor (2007) Castor, J. I. 2007, Radiation Hydrodynamics
  * Castor (2009) Castor, J. I. 2009, in Recent Directions In Astrophysical Quantitative Spectroscopy And Radiation Hydrodynamics, ed. I. Hubeny, J. M. Stone, K. MacGregor, & K. Werner (New York: American Inst. of Physics)
  * Chan & Lingenfelter (1993) Chan, K.-W., & Lingenfelter, R. E. 1993, ApJ, 405, 614
  * Chandrasekhar (1931) Chandrasekhar, S. 1931, Astrophysical Journal, 74, 81
  * Chandrasekhar (1956) —. 1956, Proceedings of the National Academy of Science, 42, 273
  * Chandrasekhar (1961) —. 1961, Hydrodynamic and hydromagnetic stability
  * Chandrasekhar & Prendergast (1956) Chandrasekhar, S., & Prendergast, K. H. 1956, Proceedings of the National Academy of Science, 42, 5
  * Choudhuri (1998) Choudhuri, A. R. 1998, The physics of fluids and plasmas : an introduction for astrophysicists /
  * Colella & Woodward (1984) Colella, P., & Woodward, P. R. 1984, Journal of Computational Physics, 54, 174
  * Collins et al. (2010) Collins, D. C., Xu, H., Norman, M. L., Li, H., & Li, S. 2010, ApJS, 186, 308
  * Cyburt et al. (2010) Cyburt, R. H., Amthor, A. M., Ferguson, R., et al. 2010, ApJS, 189, 240
  * Diamond et al. (2014) Diamond, T., Hoeflich, P., & Gerardy, C. L. 2014, ArXiv e-prints
  * Diamond et al. (2015) Diamond, T. R., Hoeflich, P., & Gerardy, C. L. 2015, ApJ, 806, 107
  * Diamond et al. (2018) Diamond, T. R., Hoeflich, P., Hsiao, E. Y., et al. 2018, ApJ, 861, 119
  * Domínguez & Hoeflich (2000) Domínguez, I., & Hoeflich, P. 2000, ApJ, 528, 854
  * Dominik (1992) Dominik, C. 1992, PhD thesis, -
  * Dominik (2009) Dominik, C. 2009, in Astronomical Society of the Pacific Conference Series, Vol. 414, Cosmic Dust - Near and Far, ed. T. Henning, E. Grün, & J. Steinacker, 494
  * Dominik & Tielens (1997) Dominik, C., & Tielens, A. G. G. M. 1997, ApJ, 480, 647
  * Fesen et al. (2007) Fesen, R. A., Hoeflich, P. A., Hamilton, A. J. S., et al. 2007, ApJ, 658, 396
  * Folatelli et al. (2013) Folatelli, G., Morrell, N., Phillips, M. M., et al. 2013, ApJ, 773, 53
  * Fransson & Sonneborn (1994) Fransson, C., & Sonneborn, G. 1994, in Frontiers of Space and Ground-based Astronomy, ed. W. Wamsteker, M. S. Longair, & Y. Kondo (Dordrecht: Kluwer), in press
  * Fryxell (2001) Fryxell, B. 2001, in American Institute of Physics Conference Series, Vol. 583, American Institute of Physics Conference Series, 310–312
  * Fryxell et al. (1991) Fryxell, B., Mueller, E., & Arnett, D. 1991, ApJ, 367, 619
  * Galbany et al. (2019) Galbany, L., Ashall, C., Hoeflich, P., et al. 2019, arXiv e-prints
  * Gall et al. (2018) Gall, C., Stritzinger, M. D., Ashall, C., et al. 2018, A&A, 611, A58
  * Gamezo et al. (2005) Gamezo, V., Khokhlov, A., & Oran, E. 2005, ApJ, 623, 337
  * Gamezo et al. (2005) Gamezo, V. N., Khokhlov, A. M., & Oran, E. S. 2005, ApJ, 623, 337
  * Gamezo et al. (2003a) Gamezo, V. N., Khokhlov, A. M., Oran, E. S., Ctchelkanova, A. Y., & Rosenberg, R. O. 2003a, Science, 299, 77
  * Gamezo et al. (2003b) —. 2003b, Science, 299, 77
  * Gardiner & Stone (2005) Gardiner, T. A., & Stone, J. M. 2005, Journal of Computational Physics, 205, 509
  * Gerardy et al. (2007) Gerardy, C. L., Meikle, W. P. S., Kotak, R., et al. 2007, ApJ, 661, 995
  * Giordano et al. (1984) Giordano, M., de Angelis, U., & Forlani, A. 1984, Journal of Physics C Solid State Physics, 17, 4089
  * Glasner et al. (2018) Glasner, S. A., Livne, E., Steinberg, E., Yalinewich, A., & Truran, J. W. 2018, MNRAS, 476, 2238
  * Gould (1971) Gould, R. J. 1971, Physica, 60, 145
  * Graboske et al. (1973) Graboske, H. C., Dewitt, H. E., Grossman, A. S., & Cooper, M. S. 1973, ApJ, 181, 457
  * Graur et al. (2016) Graur, O., Zurek, D., Shara, M. M., et al. 2016, ApJ, 819, 31
  * Hashimoto et al. (1995) Hashimoto, M., Eriguchi, Y., & Muller, E. 1995, A&A, 297, 135
  * Hillebrandt et al. (2013) Hillebrandt, W., Kromer, M., Röpke, F. K., & Ruiter, A. J. 2013, Frontiers of Physics, 8, 116
  * Hillier (1990) Hillier, D. J. 1990, A&A, 231, 116
  * Hoeflich (1990) Hoeflich, P. 1990, A&A, 229, 191
  * Hoeflich (1991) —. 1991, A&A, 246, 481
  * Hoeflich (1995a) —. 1995a, ApJ, 443, 89
  * Hoeflich (1995b) —. 1995b, ApJ, 440, 821
  * Hoeflich (2002) —. 2002, arXiv:0207103
  * Hoeflich (2003a) Hoeflich, P. 2003a, in Astronomical Society of the Pacific Conference Series, Vol. 288, Stellar Atmosphere Modeling, ed. I. Hubeny, D. Mihalas, & K. Werner, 185
  * Hoeflich (2003b) Hoeflich, P. 2003b, in Astronomical Society of the Pacific Conference Series, Vol. 288, Stellar Atmosphere Modeling, ed. I. Hubeny, D. Mihalas, & K. Werner (San Franciso: ASP), 371
  * Hoeflich (2006) —. 2006, Nuclear Physics A, 777, 579
  * Hoeflich (2017) —. 2017, Explosion Physics of Thermonuclear Supernovae and Their Signatures, ed. A. W. Alsabti & P. Murdin, 1151
  * Hoeflich et al. (2013) Hoeflich, P., Dragulin, P., Mitchell, J., et al. 2013, Frontiers of Physics, 8, 144
  * Hoeflich et al. (2004) Hoeflich, P., Gerardy, C. L., Nomoto, K., et al. 2004, ApJ, 617, 1258
  * Hoeflich et al. (2017a) Hoeflich, P., Hsiao, E. Y., Ashall, C. R., & the CSP collaboration. 2017a, ApJ, 846, 58
  * Hoeflich & Khokhlov (1996) Hoeflich, P., & Khokhlov, A. 1996, ApJ, 457, 500
  * Hoeflich et al. (1992) Hoeflich, P., Khokhlov, A., & Müller, E. 1992, A&A, 259, 549
  * Hoeflich et al. (1993) Hoeflich, P., Mueller, E., & Khokhlov, A. 1993, A&A, 268, 570
  * Hoeflich & Stein (2002) Hoeflich, P., & Stein, J. 2002, ApJ, 568, 779
  * Hoeflich et al. (1998) Hoeflich, P., Wheeler, J. C., & Thielemann, F.-K. 1998, ApJ, 495, 617
  * Hoeflich et al. (2017b) Hoeflich, P., Hsiao, E. Y., Ashall, C., et al. 2017b, ApJ, 846, 58
  * Hoeflich et al. (2017c) Hoeflich, P., Chakraborty, S., Comaskey, W., et al. 2017c, Mem. Soc. Astron. Italiana, 88, 302
  * Höflich (2002) Höflich, P. 2002, New Astronomy Reviews, 46, 475
  * Höflich (2003a) Höflich, P. 2003a, in Astronomical Society of the Pacific Conference Series, Vol. 288, Stellar Atmosphere Modeling, ed. I. Hubeny, D. Mihalas, & K. Werner, 185
  * Höflich (2003b) Höflich, P. 2003b, in Astronomical Society of the Pacific Conference Series, Vol. 288, Stellar Atmosphere Modeling, ed. I. Hubeny, D. Mihalas, & K. Werner, 371
  * Höflich (2009) Höflich, P. 2009, in American Institute of Physics Conference Series, Vol. 1171, American Institute of Physics Conference Series, ed. I. Hubeny, J. M. Stone, K. MacGregor, & K. Werner, 161–172
  * Höflich et al. (2013) Höflich, P., Dragulin, P., Mitchell, J., et al. 2013, Frontiers of Physics, 8, 144
  * Höflich et al. (2006) Höflich, P., Gerardy, C. L., Marion, H., & Quimby, R. 2006, New Astronomy Reviews, 50, 470
  * Höflich et al. (2004) Höflich, P., Gerardy, C. L., Nomoto, K., et al. 2004, ApJ, 617, 1258
  * Hoof (2010) Hoof, P. V. 2010, Atomic Line List v2.04, URL: http://www.pa.uky.edu/ peter/atomic/
  * Hoyle & Fowler (1960) Hoyle, F., & Fowler, W. A. 1960, ApJ, 132, 565
  * Hristov et al. (2018) Hristov, B., Collins, D. C., Hoeflich, P., Weatherford, C. A., & Diamond, T. R. 2018, ApJ, 858, 13
  * Hsiao et al. (2020) Hsiao, E. Y., Hoeflich, P., Ashall, C., et al. 2020, ApJ, 900, 140
  * Hubeny & Lanz (1992) Hubeny, I., & Lanz, T. 1992, A&A, 262, 501
  * Hubeny & Mihalas (2014) Hubeny, I., & Mihalas, D. 2014, Theory of Stellar Atmospheres
  * Hui (1978) Hui, A. 1978, J. Quant. Spec. Radiat. Transf., 19, 509
  * Hummer & Rybicki (1992) Hummer, D. G., & Rybicki, G. B. 1992, ApJ, 387, 248
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * Iben & Tutukov (1984) Iben, Jr., I., & Tutukov, A. V. 1984, ApJs, 54, 335
  * Isern et al. (2011) Isern, J., Hernanz, M., & José, J. 2011, in Lecture Notes in Physics, Berlin Springer Verlag, Vol. 812, Lecture Notes in Physics, Berlin Springer Verlag, ed. R. Diehl, D. H. Hartmann, & N. Prantzos, 233–308
  * Itoh et al. (1979) Itoh, N., Totsuji, H., Ichimaru, S., & Dewitt, H. E. 1979, ApJ, 234, 1079
  * Karp et al. (1977) Karp, A. H., Lasher, G., Chan, K. L., & Salpeter, E. E. 1977, ApJ, 214, 161
  * Kashi & Soker (2011) Kashi, A., & Soker, N. 2011, MNRAS, 417, 1466
  * Kerzendorf et al. (2014) Kerzendorf, W. E., Taubenberger, S., Seitenzahl, I. R., & Ruiter, A. J. 2014, ApJ, 796, L26
  * Khokhlov (2001) Khokhlov, A. 2001, APS Meeting Abstracts, C3003
  * Khokhlov (1989) Khokhlov, A. M. 1989, MNRAS, 239, 785
  * Khokhlov (1995a) —. 1995a, ApJ, 449, 695
  * Khokhlov (1995b) —. 1995b, ApJ, 449, 695
  * Khokhlov (2000) —. 2000, arXiv:0008463
  * Khokhlov et al. (1997) Khokhlov, A. M., Oran, E. S., & Wheeler, J. C. 1997, ApJ, 478, 678
  * Kromer et al. (2010) Kromer, M., Sim, S. A., Fink, M., et al. 2010, ApJ, 719, 1067
  * Kurucz (1993) Kurucz, R. 1993, CDROM No. 1: Atomic Data for Opacity Calculations, SAO Cambridge, MA
  * Kurucz & Bell (1995) Kurucz, R., & Bell, B. 1995, Atomic Line Data (R.L. Kurucz and B. Bell) Kurucz CD-ROM No. 23. Cambridge, Mass.: Smithsonian Astrophysical Observatory, 1995., 23
  * Lang (1999) Lang, K. R. 1999, Astrophysical Formulae (Berlin: Springer)
  * Langanke (2004) Langanke, K. 2004, Nuclear Astrophysics: Selected Topics, ed. J. M. Arias & M. Lozano, Vol. 652, 173
  * Leloudas et al. (2009) Leloudas, G., Stritzinger, M. D., Sollerman, J., et al. 2009, A&A, 505, 265
  * Lepp et al. (1990) Lepp, S., Dalgarno, A., & McCray, R. 1990, ApJ, 358, 262
  * Li et al. (2008) Li, S., Li, H., & Cen, R. 2008, ApJS, 174, 1
  * Liebert et al. (2003) Liebert, J., Bergeron, P., & Holberg, J. 2003, aj, 125, 348
  * Livne (1990) Livne, E. 1990, ApJ, 354, L53
  * Lorén-Aguilar et al. (2009) Lorén-Aguilar, P., Isern, J., & García-Berro, E. 2009, A&A, 500, 1193
  * Maeda et al. (2010) Maeda, K., Benetti, S., Stritzinger, M., et al. 2010, Nature, 466, 82
  * Maguire et al. (2018) Maguire, K., Sim, S. A., Shingles, L., et al. 2018, MNRAS, 477, 3567
  * Mestel (1956) Mestel, L. 1956, MNRAS, 116, 324
  * Mignone (2007) Mignone, A. 2007, J. Comput. Phys, 225, 1427
  * Mihalas et al. (1975) Mihalas, D., Kunasz, P., & Hummer, D. 1975, ApJ, 202, 465
  * Mihalas et al. (1976a) —. 1976a, ApJ, 203, 647
  * Mihalas et al. (1976b) —. 1976b, ApJ, 206, 515
  * Mihalas & Weibel Mihalas (1984) Mihalas, D., & Weibel Mihalas, B. 1984, Foundations of radiation hydrodynamics (New York: Oxford University Press)
  * Milne et al. (2001) Milne, P. A., The, L. S., & Leising, M. D. 2001, ApJ, 559, 1019
  * Motohara et al. (2006) Motohara, K., Maeda, K., Gerardy, C. L., et al. 2006, ApJl, 652, L101
  * Niemeyer et al. (1996) Niemeyer, J. C., Hillebrandt, W., & Woosley, S. E. 1996, ApJ, 471, 903
  * Niemeyer & Kerstein (1997) Niemeyer, J. C., & Kerstein, A. R. 1997, New Astronomy, 2, 239
  * Nomoto (1982) Nomoto, K. 1982, ApJ, 257, 780
  * Nomoto K. (1984) Nomoto K., Thielemann F. K., Y. K. 1984, ApJ, 286, 644
  * P.A. Milne (1999) P.A. Milne, M.D. Leising, T. 1999, ApJs, 124
  * Pakmor et al. (2012a) Pakmor, R., Edelmann, P., Röpke, F. K., & Hillebrandt, W. 2012a, MNRAS, 424, 2222
  * Pakmor et al. (2012b) —. 2012b, MNRAS, 424, 2222
  * Pakmor et al. (2011a) Pakmor, R., Hachinger, S., Röpke, F. K., & Hillebrandt, W. 2011a, A&A, 528, A117
  * Pakmor et al. (2011b) —. 2011b, A&A, 528, A117
  * Pakmor et al. (2021) Pakmor, R., Zenati, Y., Perets, H. B., & Toonen, S. 2021, MNRAS, 503, 4734
  * Patat et al. (2012) Patat, F., Hoeflich, P., Baade, D., & et.al. 2012, A&A, 545, A7
  * Penney & Hoeflich (2014) Penney, R., & Hoeflich, P. 2014, ApJ, 795, 84
  * Penny et al. (2012) Penny, B., Hoeflich, P., & Gerardy, C. 2012, ApJ
  * Petuchowski et al. (1989) Petuchowski, S. J., Dwek, E., Allen, Jr., J. E., & Nuth, III, J. A. 1989, ApJ, 342, 406
  * Phillips (1993) Phillips, M. M. 1993, Astrophysical Journal, Letters, 413, L105
  * Piersanti et al. (2003) Piersanti, L., Gagliardi, S., Iben, Jr., I., & Tornambé, A. 2003, ApJ, 598, 1229
  * Polin et al. (2019) Polin, A., Nugent, P., & Kasen, D. 2019, ApJ, 873, 84
  * Poludnenko et al. (2019) Poludnenko, A. Y., Chambers, J., Ahmed, K., Gamezo, V. N., & Taylor, B. D. 2019, Science, 366, aau7365
  * Quimby et al. (2006) Quimby, R., Hoeflich, P., Kannappan, S. J., et al. 2006, ApJ, 636, 400
  * Rappaport et al. (1994) Rappaport, S., Di Stefano, R., & Smith, J. D. 1994, Astrophysical Journal, 426, 692
  * Rasio & Shapiro (1994) Rasio, F. A., & Shapiro, S. L. 1994, ApJ, 432, 242
  * Remming & Khokhlov (2014) Remming, I. S., & Khokhlov, A. M. 2014, ApJ, 794, 87
  * Rho et al. (2020) Rho, J., Evans, A., Geballe, T. R., et al. 2020, arXiv:2010.00662
  * Röpke (2005) Röpke, F. K. 2005, A&A, 432, 969
  * Röpke et al. (2012) Röpke, F. K., Kromer, M., Seitenzahl, I. R., et al. 2012, ApJl, 750, L19
  * Scharmer (1984) Scharmer, G. B. 1984, in Methods in Radiative Transfer, ed. W. Kalkofen (Cambridge: Cambridge Univ. Press), 173
  * Schmidt et al. (2003) Schmidt, G. D., Harris, H. C., Liebert, J., et al. 2003, ApJ, 595, 1101
  * Seaton (2005) Seaton, M. J. 2005, MNRAS, 362, L1
  * Segretain et al. (1997) Segretain, L., Chabrier, G., & Mochkovitch, R. 1997, ApJ, 481, 355
  * Serge (1977) Serge, E. 1977, Nuclei and Particles (New York, NY: Benjamin and Cummings)
  * Sharp (1988) Sharp, C. M. 1988, A&AS, 72, 355
  * Sharpton et al. (1970) Sharpton, F. A., St. John, R. M., Lin, C. C., & Fajen, F. E. 1970, Phys Rev., A2, 1305
  * Shen et al. (2018) Shen, K. J., Kasen, D., Miles, B. J., & Townsley, D. M. 2018, ApJ, 854, 52
  * Shen et al. (2010) Shen, K. J., Kasen, D., Weinberg, N. N., Bildsten, L., & Scannapieco, E. 2010, ApJ, 715, 767
  * Shen & Moore (2014) Shen, K. J., & Moore, K. 2014, ApJ, 797, 46
  * Silvestri et al. (2007) Silvestri, N. M., Lemagie, M. P., Hawley, S. L., et al. 2007, AJ, 134, 741
  * Sim et al. (2010) Sim, S. A., Röpke, F. K., Hillebrandt, W., et al. 2010, ApJ, 714, L52
  * Slattery et al. (1982) Slattery, W. L., Doolen, G. D., & Dewitt, H. E. 1982, Phys. Rev. A, 26, 2255
  * Sobolev (1957) Sobolev, V. V. 1957, AZh, 34, 694
  * Sollerman et al. (2004) Sollerman, J., Lindahl, J., Kozma, C., et al. 2004, A&A, 428, 555
  * Stone et al. (1992a) Stone, J. M., Mihalas, D., & Norman, M. L. 1992a, ApJS, 80, 819
  * Stone et al. (1992b) —. 1992b, ApJS, 80, 819
  * Stritzinger et al. (2002) Stritzinger, M., Hamuy, M., Suntzeff, N. B., et al. 2002, Astronomical Journal, 124, 2100
  * Tanikawa et al. (2018) Tanikawa, A., Nomoto, K., & Nakasato, N. 2018, ApJ, 868, 90
  * Telesco et al. (2015) Telesco, C. M., Hoeflich, P., Li, D., & et.al. 2015, ApJ, 798, 93
  * Thielemann et al. (1994a) Thielemann, F., Kratz, K., Pfeiffer, B., et al. 1994a, Nuclear Physics A, 570, 329
  * Thielemann et al. (1994b) Thielemann, F., Nomoto, K., & Hashimoto, M. 1994b, in Supernovae, ed. S. A. Bludman, R. Mochkovitch, & J. Zinn-Justin, 629–+
  * Tout et al. (2008) Tout, C. A., Wickramasinghe, D. T., Liebert, J., Ferrario, L., & Pringle, J. E. 2008, MNRAS, 387, 897
  * Towns et al. (2014) Towns, J., Cockerill, T., Dahan, M., et al. 2014, Computing in Science and Engineering, 16, 62
  * Townsley et al. (2019) Townsley, D. M., Miles, B. J., Shen, K. J., & Kasen, D. 2019, arXiv:1903.10960
  * Turk et al. (2011) Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9
  * Van Der Walt et al. (2011) Van Der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, arXiv:1102.1523
  * Van Horn (1969) Van Horn, H. M. 1969, Physics Letters A, 28, 706
  * Wang et al. (2009a) Wang, B., Chen, X., Meng, X., & Han, Z. 2009a, ApJ, 701, 1540
  * Wang et al. (2009b) Wang, B., Meng, X., Chen, X., & Han, Z. 2009b, MNRAS, 395, 847
  * Wang et al. (2013) Wang, X., Wang, L., Filippenko, A. V., Zhang, T., & Zhao, X. 2013, Science, 340, 170
  * Webbink (1984) Webbink, R. F. 1984, ApJ, 277, 355
  * Whelan & Iben (1973) Whelan, J., & Iben, I. J. 1973, ApJ, 186, 1007
  * Woosley & Kasen (2011) Woosley, S. E., & Kasen, D. 2011, ApJ, 734, 38
  * Woosley & Weaver (1994) Woosley, S. E., & Weaver, T. A. 1994, ApJ, 423, 371
  * Woosley et al. (1980) Woosley, S. E., Weaver, T. A., & Taam, R. E. 1980, in Texas Workshop on Type I Supernovae, ed. J. C. Wheeler, 96–112
  * Yang et al. (2018) Yang, Y., Wang, L., Baade, D., et al. 2018, ApJ, 852, 89
  * Yang et al. (2019) Yang, Y., Hoeflich, P. A., Baade, D., et al. 2019, arXiv e-prints
  * Yoon et al. (2007) Yoon, S.-C., Podsiadlowski, P., & Rosswog, S. 2007, MNRAS, 380, 933
  * Zingale et al. (2011) Zingale, M., Nonaka, A., Almgren, A. S., et al. 2011, ApJ, 740, 8

##  Appendix A Mass and molar fractions and burning operator

Here we define the mass and molar fractions in terms of other quantities and
provide some identities. Firstly note that in a $\rm{fuel}$/$\rm{product}$
mixture, the partial mass densities add up to the total mass density:

$\rho=\rho_{\rm{fuel}}+\rho_{\rm{prod}}$ (A1)

Let the respective molar masses are $\mathcal{A}_{\rm{fuel}}$ and
$\mathcal{A}_{\rm{prod}}$. The mass, molar, and burned fractions, $X$, $Y$ and
$f$, are defined as follows:

$\displaystyle
X_{\rm{fuel}}\equiv\rho_{\rm{fuel}}/\rho,\;\;\;\;Y_{\rm{fuel}}\equiv
X_{\rm{fuel}}/\mathcal{A}_{\rm{fuel}}$ (A2) $\displaystyle
X_{\rm{prod}}\equiv\rho_{\rm{prod}}/\rho,\;\;\;\;\ Y_{\rm{prod}}\equiv
X_{\rm{prod}}/\mathcal{A}_{\rm{prod}}$ (A3) $\displaystyle f\equiv
f_{\rm{prod}}\equiv\frac{Y_{\rm{prod}}}{Y_{\rm{fuel}}+Y_{\rm{prod}}}$ (A4)
$\displaystyle f_{\rm{fuel}}\equiv 1-f$ (A5)

For the burned fraction $0\leq f\leq 1$ always holds, whereas $f=0$
corresponds to pure fuel and $f=1$ to pure product. The two mass fractions
obey $X_{\rm{fuel}}+X_{\rm{prod}}\equiv 1$.

Equation (A1) allows for a code to keep track of one partial density variable
in addition to the total mass density. In our case we had chosen the product
density. The burning operator, eq. (8), is evaluated at the end of the time
cycle. Below are the identities necessary to go back and forth between the
relevant quantities:

$\displaystyle
f=\frac{\mathcal{A}_{\rm{fuel}}\;\rho_{\rm{prod}}}{\mathcal{A}_{\rm{prod}}\;\rho+\left(\mathcal{A}_{\rm{fuel}}-\mathcal{A}_{\rm{prod}}\right)\rho_{\rm{prod}}}$
(A6)
$\displaystyle\rho_{\rm{prod}}=\frac{\mathcal{A}_{\rm{prod}}\;f}{\mathcal{A}_{\rm{fuel}}(1-f)+\mathcal{A}_{\rm{prod}}\;f}\;\rho$
(A7)

All our MHD simulations assume the product being 56Ni, as well as 12C : 16O
molar ratio of 50 : 50 for the fuel mixture. This corresponds to
$\mathcal{A}_{\rm{fuel}}:\mathcal{A}_{\rm{prod}}=1:4$. Note that $\rho_{prod}$
is a tracer field not entering the MHD-equations.

##  Appendix B Radiation Hydrodynamics with HYDRA

Our HYDrodynamical RAdiation HYDRA code consists of physics-based modules to
provide a solution for the nuclear networks, the statistical equations needed
to determine the atomic level population, the equations of state, the
opacities, and the hydro and radiation problems. The individual modules are
coupled explicitly (Fig. 14). Consistency between the solutions is achieved
iteratively by perturbation methods, a combination of accelerate lambda
acceleration plus equivalent-two level approach (ALI2) net-cooling and heating
rates for radiative coupling terms (Höflich, 2003a, 2009), and excitation and
ionization by hard radiation and particles. Different modules are employed
during different stages of the simulation. Below, we will give the basic
equations with references to the methods actually used in this paper, using
standard notation as in Mihalas & Weibel Mihalas (1984) wherever possible.
Though several of the modules are formally used as part of HYDRA but do not
affect the results in this paper as discussed in the main text.

Figure 14:  Block diagram of our numerical scheme to solve radiation
hydrodynamical problems, including detailed equation of state, nuclear and
atomic networks.

Hydrodynamics: The structure of the expanding envelopes are obtained by three
different modules by a) assuming free expansion, or by solving the non-
relativistic hydro equations in b) the Lagrangian frame for spherical geometry
including a front tracking scheme to resolve shock fronts (Fryxell, 2001), or
c) the Eulerian scheme for full 3-D using cartesian coordinates based on
PROMETHEUS (Fryxell et al., 1991). In this paper, we use module b) for the
early evolution till $\approx day10$, and module a) afterwards. The hydro
modules use an explicit Piecewise Parabolic Method (PPM) by Colella & Woodward
(1984) to solve the compressible reactive flow equations with variable
adiabatic gradients based on low and high density equation of state (see
below). PPM is implemented as a step followed by separate remaps of the
thermal and kinetic energy to avoid numerical generation of spurious pressure
disturbances during propagation of reaction fronts (flames and detonations).

Deflagration fronts: For spherical geometry, the deflagration speed in mass
coordinates is given by

$v_{t}=max(v_{cond},C_{1}\sqrt{\alpha_{T}gL_{f}})$ (B1)

with $v_{cond}$ being the conduction speed, with the Atwood number
$\alpha_{T}=(\alpha-1)/(\alpha+1)$ and
$\alpha=\rho^{+}(r_{burn})/\rho^{-}(r_{burn})$, $r_{burn}$ being the distance
of the front from the center, $\rho^{+},\rho^{-}$ the density jump across the
front. $L_{f}$ is the characteristic length scale for the freeze out of the
turbulence. The main effect of the expansion is the freeze out of the
turbulence on scales $L_{f}$ where the turbulent velocity $v_{t}$ due to RT
instabilities is comparable to the differential expansion velocities. For
these flows, $v_{t}\approx v_{exp}=L_{f}~{}t$ with $v_{exp}$ is the expansion
velocity at time $t$. $C_{1}$ has been calibrated (Domínguez & Hoeflich,
2000), based on 3D hydrodynamical simulations (Khokhlov et al., 1997; Gamezo
et al., 2003b). Here, we did use $C1=0.2$ . In other works where we employ 3D
hydro, we use the burning operator given in Equations 8 & 9.

High-density equation of state and nuclear reactions: The EOS is for a
partially degenerate and partially relativistic Bose- and Fermi-gas plus
interactions due to effects of Coulomb corrections, quantum-relativistic
effects on the electron component, and electron-positron pair production
(Nomoto K., 1984; Chandrasekhar, 1931; Van Horn, 1969; Slattery et al., 1982;
Giordano et al., 1984; Chandrasekhar, 1931; Hoeflich, 2006). Electron
screenings are taken from Graboske et al. (1973) and in the weak, intermediate
and intermediate-strong regime, and from Itoh et al. (1979) in the strong
regime.

Nuclear reactions are used in the high density regime and temperatures above
$5\times 10^{5}K$. A network is used with rates including weak, strong and
electromagnetic reactions. It is based on the implementation by Thielemann et
al. (1994a, b) but with modified matrix solvers. For isotope $i$, the change
of the abundance ratio per nucleon $Y_{i}$ is given by

$\dot{Y}_{i}=\sum_{j}N^{i}_{j}\lambda_{j}Y_{j}+\sum_{j,k}{N^{i}_{j,k}\over
1+\delta_{jk}}\rho N_{A}\langle\sigma v\rangle_{j;k}Y_{j}Y_{k}+O(3).$ (B2)

The first term on the right-hand side includes single-particle processes;
decays, photo-disintegrations, electron and positron captures, and neutrino
induced reactions with $\lambda_{j}$ being the rate per particle $j$. The
second term includes two particle reactions for particle densities
$n_{j},n_{k}$, and mass density, $\rho$, with $N_{A}$ being Avogadro’s number
and $\delta_{jk}$ is the Kronecker delta. $\langle\sigma v\rangle_{j;k}$ is
the nuclear cross-section convolved by the velocity of the particles $j$ and
$k$ in the plasma. The $N^{i}_{j}$ and $N^{i}_{j,k}$ are the number of
particles $j,k$ created or destroyed in the process. $O(3)$ includes terms due
to multi-particle reactions. Here, we use 218 isotopes, but the solver has
been used for up to 2438 isotopes. For temperatures larger than $6.5\times
10^{9}K$, we assume nuclear statistical equilibrium mediated by strong and
photon induced reactions, with only terms due to weak reactions remaining in
the rate equations. For $\langle\sigma v\rangle_{i;k}$, Boltzman distribution
is assumed for the velocity $v$ of the isotopes in the plasma. The individual
rates between isotopes are fitted by expansions in $\rho$ and $T$ with factors
tabulated in the reaction-network library REACLIB and modified weak reactions.
Updated cross-sections are as published in Cyburt et al. (2010). For more
details, see e.g. Thielemann et al. (1994a); Langanke (2004).

Low density equation of state, opacities and source functions: The data for
the atomic levels and line transitions are taken from the compilation of
Kurucz (Kurucz & Bell, 1995), Seaton (2005), and Hoof (2010), supplemented by
additional forbidden lines from Diamond et al. (2015). For the atomic models,
we reduce the number of energy levels by use of level-merging/super-levels,
with the assumption that coupling between the merged levels being in
thermodynamical equilibrium. This implies full frequency redistribution of
individual transitions (see below). For the radiation transport, Voigt
functions (Hui, 1978) are assumed for the individual line profiles in comoving
frame. For Monte Carlo transport, we assume the narrow line limit (Sobolev,
1957; Adams et al., 1971; Castor, 2004) and take into account the probability
for absorption along the way by lines of higher frequency, similar to the
formulation by Karp et al. (1977) and the narrow line-limit (Sobolev, 1957;
Hummer & Rybicki, 1992; Höflich, 2003a). The module for molecular kinetics use
rates from Petuchowski et al. (1989); Sharp (1988); Lepp et al. (1990);
Sharpton et al. (1970); Hashimoto et al. (1995); Gerardy et al. (2007); Rho et
al. (2020). In addition, the time-dependent formation of dust for carbonates,
silicates and iron-crystals has been implemented as a kinetic gas theory based
on the codes for dust formation (Dominik, 1992; Dominik & Tielens, 1997;
Dominik, 2009).

We solve the full set of statistical equations to determine the population
density of atomic states, $n_{i}$, where $i$ is shorthand for the excitation
state $i$ out of $js$ excitation states for the ionization state $k$ of
element $el$. The time dependent level populations are given by the $el\times
k\times j$ equations,

${\partial n_{i}\over\partial t}+\nabla(n_{i}{\bf v})=\sum_{i\neq
j}(n_{j}P_{ji}-n_{i}P_{ij})+n_{k^{\prime}}(R_{k^{\prime}i}+C_{k^{\prime}i})-n_{i}(R_{ik^{\prime}}+C_{ik^{\prime}}).$
(B3)

Here the rate, $P_{ij}$, is the sum of radiative, $R_{ij}$, and collisional,
$C_{ij}$, processes; $k^{\prime}$ stands for all transitions from all bound
levels of higher ionization states.

The non-thermal excitation by hard $\gamma$-rays and electrons are added into
the radiative rates $R$ with fractions according to the individual level
densities. Time scales in $\delta n_{i}/\delta t$ are dominated by the slow
bound-free transitions and, therefore, their time-dependence is taken into
account only. Thus, the time-dependent solution can be obtained from the
stationary solution $\tilde{n_{i}}$ by solving an inhomogeneous ODE
analytically or as a simple system of linear equations (Hoeflich, 1995a) of
the form

${dn_{i}\over dt}=-C_{1}n_{i}+C_{2}~{}{\rm
with}~{}C_{1}={\tilde{n}_{k}\over\tilde{n}_{i}}(R_{ki}+C_{ki}),~{}~{}{\rm
and}~{}~{}C_{2}=n_{k}(R_{ki}+C_{ki})$ (B4)

. The radiative rates between a lower and upper level $i$ and $j$,
respectively, are given by

$\tilde{R_{ij}}=4\pi\int{\alpha_{ij}(\nu)\over h\nu}J_{\nu}d\nu+R(\text{non-
thermal}),{\rm and}~{}\tilde{R_{ji}}=4\pi\int{\alpha_{ij}(\nu)\over
h\nu}\biggl{[}{{2h\nu^{3}\over c^{2}}+J_{\nu}}e^{-{h\nu\over
kT}}\biggl{]}d\nu.$ (B5)

with $R_{ij}=(n_{j}/n_{i})^{*}\tilde{R_{ij}}$ and $R(non-thermal)$ being the
energy input by hard-radiation and non-thermal electrons, $J_{\nu}$ is the
first momentum of the radiation field, $\alpha_{ij}$ is the cross-section, and
$*$ denotes values for thermodynamical equilibrium.

The equation of particle and charge conservation can be written in the
following form

$\sum_{el}\sum_{k}\biggl{[}\sum_{levels}n_{el,k,j}\biggr{]}=n_{tot}{\rm~{}~{}and~{}~{}}\sum_{el}\sum_{k}Z_{k}\biggl{[}\sum_{levels}n_{el,k,j}\biggr{]}=n_{e}$
(B6)

where $el$, $k$ and $level$ are the sums of the element $el$, ionization
stages $k$ for element and levels in the particular ion with Zk being the
charge of the ion. The electron and total densities are given by $n_{e}$ and
$n_{tot}$, respectively. Complete redistribution over each individual line
both in frequency and in angle is assumed in the comoving frame. Complete
redistribution also implies that the relative populations within the sub-
levels or the merged levels are described by a Maxwell-Boltzmann distribution.

The formal total source function $S_{\nu}$, emissivity $\eta$, the opacity
$\chi$ and frequency redistribution function $\psi$ are related by

$S_{\nu}={\eta_{\nu}\over\chi_{\nu}}\psi_{\nu}$ (B7)

with

$\eta_{\nu}={2h\nu^{3}\over
c^{2}}\sum_{el}\sum_{k}\Biggl{[}\sum^{I}_{i=1}\biggl{[}n^{*}_{i}\alpha_{ik}(\nu)e^{-h\nu\over
kT}+\sum^{I}_{j>i}(\alpha_{ij}(\nu)n_{j}{g_{i}\over
g_{j}})\biggr{]}+n_{e}n_{ion}\alpha_{ff}(\nu)e^{-h\nu/kT}\Biggr{]}+\Biggl{[}n_{e}\sigma_{e}+\sum_{\ell^{*}}\chi_{\ell^{*}}\Biggr{]}J_{\nu}$
(B8)

and

$\chi_{\nu}=\sum_{el}\sum_{k}\Biggl{[}\sum^{I}_{i=1}\Bigl{[}\alpha_{ik}(\nu)(n_{i}-n^{*}_{i}e^{-h\nu\over
kT})+\sum^{I}_{j>i}\alpha_{ij}(\nu)\Bigl{(}n_{i}-n_{j}{g_{i}\over
g_{j}}\Bigr{)}\biggr{]}+n_{e}n_{ion}\alpha_{ff}(\nu)(1-e^{-h\nu\over
kT})\Biggr{]}+\Biggl{[}n_{e}\sigma_{e}+\sum_{\ell^{*}}\chi_{\ell^{*}}\Biggr{]}$
(B9)

with $\alpha$ being the cross sections for bound-bound, bound-free and free-
free transitions. $\chi_{\ell^{*}}$ are opacities of those weak scattering
lines not treated in non-LTE.

Line redistribution functions: For the source function of low energy photons,
complete redistribution is assumed in general for each individual bound-bound
transition. However, due to the large number of lines, they will overlap even
in comoving frame due to their natural line width which results in a frequency
redistribution in case of radiation transport in spherical geometry 222as used
in LC simulations. The redistribution function $\psi_{\nu}$ is evaluated
numerically by weighting the neighboring frequencies with the overlap of a
specific transition so that the frequency integral over $S_{\nu}$ is conserved
(Adams et al., 1971; Hummer & Rybicki, 1992; Rappaport et al., 1994; Hoeflich,
1995a).333Note that the frequency coupling stabilize the iteration scheme.

Coupling of radiation transport, statistical and hydro equations: We use the
well established method of accelerated lambda iteration (ALI, e.g. Cannon
(1973); Scharmer (1984); Hillier (1990); Hoeflich (1990); Hubeny & Lanz
(1992)). We employ several concepts to improve the stability, and convergence
rate/control including the concept of leading elements, the use of net rates,
level locking, reconstruction of global photon redistribution functions,
equivalent-2-level approach (Athay, 1972; Avrett & Loeser, 1988; Hoeflich,
1995a), and predictor-corrector methods.

Radiation transport for low energy photons: For the time-dependent transport,
we use variable Eddington Tensor solvers, which are implicit in time (Mihalas
& Weibel Mihalas, 1984; Stone et al., 1992a; Höflich, 2003a; Castor, 2007,
2009; Höflich, 2009)444Some new additional VET modules are in the verification
phase, and stationary solutions of the transport for the closures.

For spherical geometry, we solve the following equations

${1\over r^{2}}{\partial(r^{2}H_{\nu})\over\partial r}-{\nu_{o}v(r)\over
cr}\biggl{[}{\partial(1-f_{\nu})J_{\nu}\over\partial\nu}+\beta(r){\partial
f_{\nu}J_{\nu}\over\partial\nu}\biggr{]}=\eta_{\nu}-\chi_{\nu}J_{\nu}$ (B10)

and

${\partial f_{\nu}J_{\nu}\over\partial r}+{3f_{\nu}-1\over
r}J_{\nu}-{\nu_{o}v(r)\over
cr}\biggl{[}{\partial(1-g_{\nu})H_{\nu}\over\partial\nu}+\beta(r){\partial(g_{\nu}H_{\nu})\over\partial\nu}\biggr{]}=-\chi_{\nu}H_{\nu},$
(B11)

where $\beta$ and the Eddington factors are defined in the usual way

$\beta={d~{}ln(v(r))\over d~{}ln(r)},~{}~{}f_{\nu}=K_{\nu}/J_{\nu},{\rm
and}~{}g_{\nu}=N_{\nu}/H_{\nu}.$ (B12)

$J_{\nu},H_{\nu},K_{\nu}$ and $N_{\nu}$ are the first four moments of the
intensity. The 4th moment is needed as a closure relation of the system. The
Eddington factors $f_{\nu}$ and $g_{\nu}$ are obtained from solutions for the
stationary case by integration along rays in a Rybicki-like scheme for non-
relativistic and relativistic velocity fields including advection terms
((Mihalas et al., 1975, 1976a, 1976b), MKH methods). Note that the Eddington
factors includes the frequency redistribution function $\psi$. Time
independent Eddington factors are assumed during each time step.

For the 3D case, the variable Eddington-Tensor method is used which is
accurate to the order O(u/v) and we neglect acceleration terms (Buchler, 1979,
1983; Stone et al., 1992b; Höflich, 2009; Hubeny & Mihalas, 2014). The system
of equations is given by

${1\over c}{\delta F_{\nu}\over\delta t}+{1\over
c}\nabla(uF_{\nu}+c\nabla(T_{\nu}E_{\nu}))=-k_{\nu}F_{\nu}$ (B13)

$\rho{D(E_{\nu}/\rho)\over Dt}+E_{\nu}T_{\nu}\nabla
u-\nabla[{cT_{\nu}\over\kappa_{n}u\rho q_{\nu}}\nabla(q_{\nu}E_{\nu})]=4\pi
j_{\nu}-\kappa_{\nu}\rho cE_{\nu}$ (B14)

with the tensor T and, following Auer (1971), the variable $q_{\nu}$ as
defined by:

$\displaystyle\nabla{T_{\nu}\nabla log~{}q_{\nu}}=\nabla(\nabla T_{\nu})$
(B15) $\displaystyle{\partial ln(r^{2}q_{\nu})\over\partial
r}={3f_{\nu}-1\over f_{\nu}r}$ (B16)

and

$E_{\nu}(\vec{x})=1/c\int\vec{I_{\nu}}(\vec{x},\vec{\Omega})d\vec{\Omega},~{}~{}\vec{F}_{\nu}(\vec{x})=1/c\int\vec{I_{\nu}}(\vec{x},\vec{\Omega})\vec{\Omega}\vec{d\Omega},~{}{\rm
and}~{}~{}\underline{P}_{\nu}(\vec{x})=1/c\int\vec{I_{\nu}}(\vec{x},\vec{\Omega})\vec{\Omega}\otimes\vec{\Omega}d\vec{\Omega}.$
(B17)

with the integrals over $\nu$ denote the corresponding radiative components in
the hydro-equations. For spherical geometry, the pressure is taken as the
diagonal elements of $\underline{P}$. For 3D, the advection frequency
redistribution terms enter via the Monte Carlo (MC) solution. For some tests
of 3D vs. spherical solutions, see Hoeflich (2002). Again, we use stationary
transport for the closure relation and tested the result against the spherical
transport module. MC methods are used including advection and aberration for
the closure relation of low energy photons for flux and polarization spectra,
and similar to hard $\gamma$ and positron transport (Hoeflich, 1991; Hoeflich
et al., 1992; Hoeflich, 1995b; Penney & Hoeflich, 2014).

Gamma-ray and positron transport: The $\gamma$-ray transport is computed in
multi-dimensions using a Monte Carlo method (Hoeflich et al., 1992; Hoeflich,
2003a) including relativistic effects and a consistent treatment of continuum
and line opacities. Interactions are calculated in the local rest frame by
transformation between observer and comoving frame. Photon and positron
packages may persist until the next hydro timestep. If the photon/positron
travel time exceeds the timestep, they’re fed to the next timestep. This
persists until they are scattered down to low energies, X-rays, or escape the
computational domain. Each package keeps a time-tracer that monitors its
travel time.

The interaction processes allowed are: Compton scattering according to the
full angle-dependent Klein-Nishina formula, pair-production, and bound-free
transitions (Ambwani & Sutherland, 1988), and the NIST XCOM data base (Berger
et al., 1998).

###  B.1 Positron Creation and Cross Sections:

The primary source of positrons is the $\beta^{+}$ channel of the
${}^{56}Co\rightarrow^{56}Fe^{*}$ which accounts for about $18\%$ of all
${}^{56}Co$ decays. About $1.4MeV$ of the total excitation energy of
${}^{56}Fe^{*}$ is available for the positrons with an energy spectrum (Serge,
1977) given by

$N(E)=Cp^{2}(E_{o}-E)^{2}(2\pi\eta(1-\exp(-2\pi\eta))^{-1})$ (B18)

where $C$ is a scale factor and $\eta$ is the charge of the nucleus times
$\hbar$ over the velocity of the electron. The mean energy of the spectrum is
.44 MeV.

For the Monte Carlo transport, we take into account three processes for the
interaction: scattering on electrons, scattering on nuclei in a plasma, and
annihilation $e^{+}e^{-}\rightarrow 2\gamma$ or $e^{+}e^{-}\rightarrow
3\gamma$ for para- and ortho-positronium, respectively.

The annihilation cross-section (Lang, 1999) is given by

$\sigma=\frac{\pi
r_{o}}{\gamma+1}[\frac{\gamma^{2}+4\gamma+1}{\gamma^{2}-1}ln(\gamma+\sqrt{\gamma^{2}-1})-\frac{\gamma+3}{\sqrt{\gamma^{2}-1}}]$
(B19)

with $r_{o}=e^{2}/mc^{2}$ and $\gamma$ is the Lorentz factor.

Following Chan & Lingenfelter (1993) and Gould (1971), the energy loss by
interaction with charged particles is given by

$\frac{dE}{dx}=\frac{-4\pi
r_{o}^{2}m_{e}c^{2}q}{AM_{n}\beta^{2}}(qln(\frac{\sqrt{\gamma-1}\gamma\beta
m_{e}c^{2}}{b_{max}})\Pi(\gamma)$ (B20)

$\Pi(\gamma)$ is the relativistic correction given by

$\Pi(\gamma)=\frac{\beta^{2}}{12}[\frac{23}{2}+\frac{7}{\gamma+1}+\frac{5}{(\gamma+1)^{2}}+\frac{2}{(\gamma+1)^{3}}].$
(B21)

$q$ is the relative charge of the particles in the media. For electron
interactions, $q$ is the atomic number $Z$ for atom. For plasma scattering,
$q$ is the ionization fraction. $b_{max}$ is the maximum impact parameter. It
is the maximum amount of energy the positron can lose in one interaction. For
electron scattering, is it the ionization potential. For the case of ions in
the plasma, the maximum impact parameter is $\hbar\omega$, which is the plasma
frequency as an impact with an ion sets up a disturbance in the plasma whose
energy is proportional to the frequency. Upon undergoing either type of
interaction, the particles are emitted isotropically in the comoving frame.

The positron transport (Penney & Hoeflich, 2014) is solved via a Monte Carlo
method very similar to our photon transport but the integration is along a
spiral path imposed by the local Lorentz force.

  *[SNe]: supernovae
  *[]:
  *[SNe Ia]:   () Type Ia
  *[DDT]: deflagration-to-detonation transition
  *[WD]: white dwarf
  *[WDs]: white dwarf
  *[RT]: Rayleigh-Taylor
  *[NIR]: near infrared
  *[MHD]: magnetohydrodynamics
  *[non-LTE]: non-local thermodynamic equilibrium
  *[IR]: infrared
  *[LCs]: light curve
  *[LC]: light curve
  *[MIR]: mid infrared
  *[QSE]: quasi nuclear statistical equilibrium