# Core-collapse supernovae simulations with reduced nucleosynthesis networks

Gerard Navó Institut für Kernphysik, Technische Universität Darmstadt,
Schlossgartenstr. 2, Darmstadt 64289, Germany Gerard Navó
<EMAIL_ADDRESS>Moritz Reichert Departament
d’Astronomia i Astrofísica, Universitat de València,
Edifici d´Investigació Jeroni Munyoz, C/ Dr. Moliner, 50, E-46100 Burjassot
(València), Spain  Martin Obergaulinger Departament d’Astronomia i
Astrofísica, Universitat de València,
Edifici d´Investigació Jeroni Munyoz, C/ Dr. Moliner, 50, E-46100 Burjassot
(València), Spain  Almudena Arcones Institut für Kernphysik, Technische
Universität Darmstadt, Schlossgartenstr. 2, Darmstadt 64289, Germany  GSI
Helmholtzzentrum für Schwerionenforschung GmbH, Planckstr. 1, Darmstadt 64291,
Germany

###### Abstract

We present core-collapse supernovae simulations including nuclear reaction
networks which impact explosion dynamics and nucleosynthesis. The different
composition treatment can lead to changes in the neutrino heating in the
vicinity of the shock, by modifying the amount of nucleons and thus the
$\mathrm{\nu}$-opacity of the region. This reduces the ram pressure outside
the shock and allows an easier expansion. The energy released by the nuclear
reactions during collapse also slows down the accretion, and aids the shock
expansion. In addition, nuclear energy generation in the post-shocked matter
produces more energetic explosions, up to $20\,\%$. Nucleosynthesis is
affected due to the different dynamic evolution of the explosion. Our results
indicate that the energy generation from nuclear reactions helps to sustain
late outflows from the vicinity of the proto-neutron star (PNS), synthesizing
more neutron-rich species. Furthermore, we show that there are systematic
discrepancies between the ejecta calculated with in-situ and ex-situ reaction
networks. The mass fractions of some Ca, Ti, Cr, and Fe isotopes are
consistently under-produced in post-processing calculations, leading to
different nucleosynthesis paths. Therefore, large in-situ nuclear reaction
networks are needed for a more accurate nucleosynthesis.

Core-collapse supernovae(304) — Supernovae(1668) — Explosive
nucleosynthesis(503) — Supernova dynamics(1664) — Astrophysical explosive
burning(100) — Nuclear astrophysics(1129)

††journal: ApJ

##  1 Introduction

Massive stars ($M\gtrsim 8\,\mathrm{M_{\odot}}$) usually end their life in
core-collapse supernovae (CCSN), explosive events that are triggered by the
collapse of the iron core and release several solar masses of products
synthesized during the stellar evolution and the explosion itself into the
interstellar medium. Thus, they play a critical role in the chemical history
of the universe and have been matter of study for many decades. In order to
investigate the mechanisms involved, radiation-hydrodynamic simulations have
been performed (for extended reviews see, e.g., Bethe, 1990; Kotake et al.,
2006; Janka, 2012; Müller, 2020; Burrows & Vartanyan, 2021), following
different evolutionary stages from collapse to shock breakout. Nevertheless,
the whole puzzle is still far from being completed.

In recent years, significant advances have been reported in multidimensional
simulations (e.g., Lentz et al., 2015; Janka et al., 2016; Roberts et al.,
2016; O’Connor & Couch, 2018; Obergaulinger & Aloy, 2020; Burrows et al.,
2020; Kuroda et al., 2020; Sandoval et al., 2021; Nakamura et al., 2022),
impact of magnetic fields (e.g., Mösta et al., 2015; Obergaulinger et al.,
2018; Bugli et al., 2021; Varma et al., 2022), pre-supernova models (e.g.,
Müller et al., 2017; Fields & Couch, 2021; Yoshida et al., 2021; Vartanyan et
al., 2021), neutrino transport (for an extended review, see Mezzacappa et al.,
2020) and reactions (e.g., Balasi et al., 2015, and references therein), high-
density equations of state (EOS) (e.g., Schneider et al., 2019; Yasin et al.,
2020) or nucleosynthesis calculations (e.g., Eichler et al., 2017; Wanajo et
al., 2018; Curtis et al., 2019; Sieverding et al., 2020; Witt et al., 2021;
Reichert et al., 2021). In this work, we focus on the impact of the different
ways to treat nuclear reactions and the associated energy generation on the
simulations.

A nuclear reaction network evolves the abundances of a set of $N$ species by a
system of $N$ coupled ordinary differential equations including all the
reactions between the species. Therefore, its size depends on the nature of
the environment. In CCSNe, it can include of the order of several hundreds
(e.g., Woosley & Weaver, 1995; Rauscher et al., 2002). Thus, including nuclear
reaction networks in multidimensional CCSNe simulations is very challenging
because of the computational cost of evolving them together with the
hydrodynamics. Therefore, hydrodynamic simulations usually employ a simple
treatment of the composition. A detailed nucleosynthesis can be obtained by
post-processing (e.g., Winteler et al., 2012; Lippuner & Roberts, 2017) using
Lagrangian tracer particles.

If temperatures are high enough ($T\gtrsim 5-6\,\mathrm{GK}$), fusion
reactions and photodisintegrations are in so called nuclear statistical
equilibrium (NSE). Within this equilibrium, the production rate of a specific
nuclear species equals its destruction rate and for constant thermodynamic
conditions the composition is also constant over time. In NSE, the abundances
are direct functions of the density, temperature, and electron fraction
($Y_{\mathrm{e}}$) conditions. Thus, instead of the costly evolution of the
reaction network, they can be used at high temperatures and included in high-
density EOS tables. If the temperature is lower, NSE does not hold anymore and
other methods are necessary to estimate the composition and energy generation.
In this case, other simplified treatments are used. The so-called flashing
scheme (Rampp & Janka, 2002) takes into account neutrons, protons,
$\mathrm{\alpha}$-particles, and a characteristic nucleus which depends on the
thermodynamic conditions, considering instantaneous burning. For example, when
silicon burning occurs, the silicon mass fraction ($X(^{28}\mathrm{Si})$) is
immediately transformed to ${}^{56}\mathrm{Ni}$ and set to
$X(^{28}\mathrm{Si})=0$. This instantaneous transformation of nuclei prevents
an inclusion and feedback of the nuclear energy generation to the
hydrodynamical simulations. An alternative is to use reduced reaction networks
(e.g., Müller, 1986; Hix & Thielemann, 1999; Timmes et al., 2000) in order to
track a small set of nuclei together with the hydrodynamics. The considered
nuclei are chosen to involve the main reactions that release or consume
internal energy and to represent the main contributions to the baryonic part
of the pressure and neutrino opacities (Cernohorsky & van Weert, 1992). In
practice, $\alpha$-chains are the most commonly used for this purpose.

Multiple groups have employed reduced networks in their works (e.g., Bruenn et
al., 2016, 2020; Nakamura et al., 2014; Couch et al., 2015; Harris et al.,
2017; Wongwathanarat et al., 2017; Sandoval et al., 2021). Bruenn et al.
(2006) mentioned that when including a $14$ $\mathrm{\alpha}$-chain in the
hydrodynamics, nuclear burning in the oxygen layer deposited extra pressure in
the vicinity of the shock and assisted its expansion. Nakamura et al. (2014)
performed simulations with a simplified light-bulb neutrino treatment using a
$13$ isotope $\mathrm{\alpha}$-network to study how the energy released by
nuclear reactions affect the dynamics of the explosion. They conclude that
energy produced by the infalling material behind the shock could aid the
explosion, especially in models with marginal explosions. A large network was
included for the first time by Harris et al. (2017) in 2D simulations with
accurate $\mathrm{\nu}$-transport. They used a $14$ $\mathrm{\alpha}$-chain
and a $150$ isotope network in axisymmetric (2D) models and studied the
uncertainties derived from post-processing nucleosynthesis. Their results
showed the limitations of using post-processing Lagrangian tracer particles
and support including reaction networks in the simulations.

In this work, the aforementioned studies motivated us to investigate in
detail, in a state-of-the-art CCSN code, how the different treatments of the
composition employed in the simulations at low temperatures have an impact in
the dynamics of the explosion and the nucleosynthesis outcomes. In Section 2,
we introduce the radiation-hydrodynamics code we used for the simulations
(Section 2.1), the reduced network module implementation as well as the chosen
networks (Section 2.2), and the different models (Section 2.3). In Section 3,
we discuss the results and in Section 4, we close with our conclusions and
final remarks.

##  2 Methods

###  2.1 Hydrodynamic code

The simulations are performed using the code Aenus-Alcar (Just et al., 2015;
Obergaulinger & Aloy, 2017), which combines special relativistic
magnetohydrodynamics (MHD), a pseudo relativistic gravitational potential, and
state of the art two-moment (M1) neutrino transport.

Prior to the implementation of the reduced network module, Aenus-Alcar
employed for matter at high densities ($\rho>\rho_{\mathrm{th}}\sim
10^{7-8}\mathrm{g}/\mathrm{cm}^{3}$) a nuclear EOS assuming NSE between
protons, neutrons, alpha particles, and a representative heavy nucleus. At
lower densities, ($\rho\leq\rho_{\mathrm{th}}$), a Helmholtz type EOS was
used, which takes into account leptonic, photonic, and baryonic contributions.
For the latter, a version of the ”flashing” scheme (Rampp & Janka, 2002)
calculated the composition of the neutrons, protons, and a characteristic
nucleus that, depending on the thermodynamic conditions, is
${}^{28}\mathrm{Si}$ or ${}^{56}\mathrm{Ni}$. In this work, we have added a
nuclear reaction network module to describe the composition outside of the NSE
regime more accurately.

###  2.2 Reduced Network Implementation

####  2.2.1 ReNet

ReNet is a highly flexible nuclear reaction network code that allows to
calculate abundance flows for networks of different sizes and complexities.
Similar to previous works in the literature (Müller, 1986; Hix & Thielemann,
1999; Timmes et al., 2000) it solves the set of ordinary first order
differential equations in an implicit way. The computational cost of nuclear
networks can be reduced by taking advantage of sparse matrices (e.g., Hix &
Thielemann, 1999; Lippuner & Roberts, 2017). However, for the reduced reaction
networks used here, the computational overhead for this is large and
consequently we use solvers for dense matrices. Within ReNet we furthermore
assume some species to be in steady state in order to consider them indirectly
and avoiding to use a differential equation for them. This allows us to save
computational cost (for more information see Appendix A and c.f., Weaver et
al., 1978; Timmes, 1999; Paxton et al., 2011, for a similar approach).

Additionally, we compute the energy released by the nuclear reactions through

$\dot{E}_{\text{nuc}}=-\sum_{i}N_{A}{\Delta m}_{i}c^{2}\dot{Y}_{i},$ (1)

where $\dot{E}_{\text{nuc}}$ is the variation in the internal energy due to
the nuclear reactions, ${\Delta m}_{i}$ is the mass excess of each specie, and
$\dot{Y}_{i}$ is the abundance time variation of the nucleus $i$.

####  2.2.2 ReNet in Aenus-Alcar

As mentioned above, previous versions of Aenus-Alcar distinguish between two
different density ranges where two different EOS regimes, one of which assumes
NSE, are used. In contrast, for the composition, the switch between NSE and
the non-NSE regime using the reaction network is based on temperature rather
than density. We found numerical artifacts in regions where the two criteria
contradicted each other, e.g., zones of low density, but high temperature for
which the EOS would not select NSE, but the composition would be given by NSE
instead of the network. This is the case for a large number of zones in the
top left part of the phase diagram shown for a typical CCSN simulation in
Figure 1. To avoid such an inconsistency, we use the same temperature criteria
in the NSE-network and the EOS transitions, such that above a certain
temperature threshold, $T>T_{\mathrm{th}}\sim 5-6\,\mathrm{GK}$, NSE is
assumed both for the composition and the EOS. Otherwise, the network is
employed together with the subnuclear EOS. Though in principle, the latter
case would also apply to very cold zones with densities up to and beyond the
nuclear saturation density, in practice this regime (bottom right part of the
figure) is not achieved in CCSN and thus of no concern.

Figure 1: Density and temperature achieved in a characteristic CCSN
simulation. It focuses on the region of the $\rho$-$T$ plane close to the
transition between the two EOS. The color depicts the amount of mass that
fulfilled each $(\rho,T)$ condition in the run. The dashed line indicates the
temperature threshold criterion introduced in this work. Finally, the dotted
line corresponds to the density threshold criterion present in previous
versions of Aenus-Alcar.

In order to prevent discontinuities in the thermodynamic quantities such as
the pressure or the internal energy, we interpolate them linearly between the
two EOSs in the temperature regime $T_{\mathrm{th}}>T>T_{\mathrm{inter}}$.
Furthermore, at $T=T_{\mathrm{th}}$, the initial values for the network need
to be given. Such seeds are provided by an NSE solver which considers the same
nuclei as the nuclear reaction network. Once the temperature in a grid cell
drops below $T_{\mathrm{th}}$, the network module evolves the composition
further until the end of the simulation or until the conditions for NSE are
reached again.

The energy from the nuclear reactions (Eq.(1)) is coupled as a source term to
the hydrodynamics energy equation. In order to avoid thermodynamic
inconsistencies, it is done for temperatures $T\leq T_{\mathrm{inter}}$. We
consider two different ways to include it. On the one hand, the feedback of
the nuclear energy from the reaction network is considered at every cell with
$T\leq T_{\mathrm{inter}}$. On the other hand, the nuclear energy generation
is taken into account only after the shock has passed the cell, i.e., in the
post-shocked region. For that purpose, we use a simple shock-detector
algorithm. We additionally consider the case with no feedback of the nuclear
energy generation from the nuclear reaction network.

####  2.2.3 Networks

Table 1: List of RN16 nuclei

RN16 Network
---
$\mathrm{n}$ | ${}^{1}\mathrm{H}$ | ${}^{4}\mathrm{He}$ | ${}^{12}\mathrm{C}$ | ${}^{16}\mathrm{O}$ | ${}^{20}\mathrm{Ne}$ | ${}^{24}\mathrm{Mg}$
${}^{27}\mathrm{Al}^{\dagger}$ | ${}^{28}\mathrm{Si}$ | ${}^{31}\mathrm{P}^{\dagger}$ | ${}^{32}\mathrm{S}$ | ${}^{35}\mathrm{Cl}^{\dagger}$ | ${}^{36}\mathrm{Ar}$ | ${}^{39}\mathrm{K}^{\dagger}$
${}^{40}\mathrm{Ca}$ | ${}^{43}\mathrm{Sc}^{\dagger}$ | ${}^{44}\mathrm{Ti}$ | ${}^{47}\mathrm{V}^{\dagger}$ | ${}^{48}\mathrm{Cr}$ | ${}^{51}\mathrm{Mn}^{\dagger}$ | ${}^{52}\mathrm{Fe}$
${}^{53}\mathrm{Fe}^{\dagger}$ | ${}^{54}\mathrm{Fe}$ | ${}^{55}\mathrm{Co}^{\dagger}$ | ${}^{56}\mathrm{Ni}$ |  |  |

Note. — Steady state nuclei are marked with a dagger.

Table 2: List of RN94eq148 nuclei

RN94eq148 Network
---
$\mathrm{n}$ | ${}^{1}\mathrm{H}$ | ${}^{4}\mathrm{He}$ | ${}^{12}\mathrm{C}$
${}^{16}\mathrm{O}$ | ${}^{20}\mathrm{Ne}$ | ${}^{23}\mathrm{Mg}^{\dagger}$ | ${}^{24}\mathrm{Mg}$
${}^{27}\mathrm{Al}^{\dagger}$ | ${}^{28,30}\mathrm{Si}$ | ${}^{27,29}\mathrm{Si}^{\dagger}$ | ${}^{29,31}\mathrm{P}^{\dagger}$
${}^{30-34}\mathrm{S}$ | ${}^{31,33}\mathrm{S}^{\dagger}$ | ${}^{33,35}\mathrm{Cl}^{\dagger}$ | ${}^{34-38}\mathrm{Ar}$
${}^{35,37}\mathrm{Ar}^{\dagger}$ | ${}^{37,39}\mathrm{K}^{\dagger}$ | ${}^{38-42,48-52}\mathrm{Ca}$ | ${}^{39,41,49,51}\mathrm{Ca}^{\dagger}$
${}^{41,43,49-55}\mathrm{Sc}^{\dagger}$ | ${}^{42-56}\mathrm{Ti}$ | ${}^{43-55}\mathrm{Ti}^{\dagger}$ | ${}^{45-57}\mathrm{V}^{\dagger}$
${}^{46-58}\mathrm{Cr}$ | ${}^{47-57}\mathrm{Cr}^{\dagger}$ | ${}^{49-61}\mathrm{Mn}^{\dagger}$ | ${}^{50-64}\mathrm{Fe}$
${}^{51-63}\mathrm{Fe}^{\dagger}$ | ${}^{53-67}\mathrm{Co}^{\dagger}$ | ${}^{54-70}\mathrm{Ni}$ | ${}^{55-69}\mathrm{Ni}^{\dagger}$
${}^{61-73}\mathrm{Cu}^{\dagger}$ | ${}^{60-74}\mathrm{Zn}$ | ${}^{59-73}\mathrm{Zn}^{\dagger}$ | ${}^{63-77}\mathrm{Ga}^{\dagger}$
${}^{64-80}\mathrm{Ge}$ | ${}^{65-79}\mathrm{Ge}^{\dagger}$ | ${}^{67-83}\mathrm{As}^{\dagger}$ | ${}^{68-84}\mathrm{Se}$
${}^{69-83}\mathrm{Se}^{\dagger}$ | ${}^{75-85}\mathrm{Br}^{\dagger}$ | ${}^{76-86}\mathrm{Kr}$ | ${}^{77-85}\mathrm{Kr}^{\dagger}$
${}^{79-87}\mathrm{Rb}^{\dagger}$ | ${}^{80-88}\mathrm{Sr}$ | ${}^{81-87}\mathrm{Sr}^{\dagger}$ | ${}^{85-89}\mathrm{Y}^{\dagger}$
${}^{86-90}\mathrm{Zr}$ | ${}^{87,89}\mathrm{Zr}^{\dagger}$ | ${}^{91}\mathrm{Nb}^{\dagger}$ | ${}^{92}\mathrm{Mo}$

Note. — Steady state nuclei are marked with a dagger. We use the ${}^{i-j}X$
notation to account for the species in
$\\{{{}^{i}X},{{}^{i+2}X},...,{{}^{j}X}\\}$.

In this work we implemented two reaction network configurations into Aenus-
Alcar. On the one hand, we have chosen a $16\text{-}\alpha$ chain (RN16),
similar to the widely used $19$ isotope approximation network from Weaver et
al. (1978), which approximately reproduces the nuclear energy generation
within CCSNe simulations and therefore provides feedback to the total energy
of the system. On the other hand, we have selected a 94-species network
(RN94eq148), which is able to track the main species synthesized in the CCSN
(Eichler et al., 2017; Harris et al., 2017), with $Y_{e}$ between $0.4$ and
$0.6$. Hence, it allows us to track slightly neutron and proton rich
trajectories along stability up to ${}^{92}\mathrm{Mo}$. In addition, $148$
species are considered in steady state, reproducing a network consisting of
$242$ nuclei. Therefore, we present the most complete network up to date,
which is evolved in CCSN simulations with state-of-the-art M1 neutrino
transport. The species included in both networks are listed in Table 1 and
Table 2.

For comparison, the nucleosynthesis is additionally computed in post-
processing with the full nuclear reaction network WinNet (Winteler et al.,
2012) and the reduced networks RN16 and RN94eq148. We consider the unbound
matter, i.e., matter with positive total energy and positive radial velocity,
at $t_{\mathrm{f}}=1.5\,\mathrm{s}$ of simulation. The evolution of the ejecta
is followed by Lagrangian tracer particles calculated backwards in time as
described in Reichert et al. (2022). At the final simulation time, we place
them at random positions in all cells flagged as unbound. The total mass
contained in a cell is distributed equally among its tracers. Their number is
set such that they have a maximum mass, $M=10^{-4}\,\mathrm{M_{\odot}}$, and
each cell contains at least $4$ tracers. In the end, the nucleosynthesis of
each tracer particle is weighted with its corresponding mass.

Table 3: List of the different models

Model | Comp | $\dot{E}_{\text{nuc}}$ | $Q^{\mathrm{LD}}_{\mathrm{\nu}}$ | Dim
---|---|---|---|---
1D_flsh | flashing | - | yes | 1d
1D_flsh* | flashing | - | no | 1d
2D_flsh | flashing | - | yes | 2d
1D_NSE | NSE | - | yes | 1d
1D_NSE* | NSE | - | no | 1d
1D_RN16 | RN16 | - | yes | 1d
1D_RN16* | RN16 | - | no | 1d
1D_RN16e | RN16 | a.s | yes | 1d
1D_RN16E | RN16 | yes | yes | 1d
2D_RN16E | RN16 | yes | yes | 2d
1D_RN94 | RN94eq148 | - | yes | 1d
1D_RN94* | RN94eq148 | - | no | 1d
1D_RN94e | RN94eq148 | a.s | yes | 1d
1D_RN94E | RN94eq148 | yes | yes | 1d
2D_RN94E | RN94eq148 | yes | yes | 2d

Note. — The second column shows the treatment of the composition. The third
column specifies whether the energy generation from the nuclear reactions is
taken into account. a.s indicates that is switched on only behind the shock.
The fourth column gives the density threshold above which the heating source
term of the neutrinos is considered. Finally, the last column states for the
dimensionality of the simulations.

###  2.3 Models

We have performed simulations ($t_{\mathrm{f}}=1.5\,\mathrm{s}$) using the
solar metallicity $20\,\mathrm{M_{\odot}}$ mass progenitor from Woosley &
Heger (2007). We map its pre-collapse composition to those of our two
networks. We consider the EOS transition at
$T_{\mathrm{th}}=5.8\,\mathrm{GK}$, given that NSE breaks down around this
temperature, and $T_{\mathrm{inter}}=5.0\,\mathrm{GK}$. Since this work is
focused on the impact of the composition treatment at $T<T_{\mathrm{th}}$, all
of the models have the same configuration at high temperatures. We apply the
SFHo EOS (Steiner et al., 2013) and assume the composition in NSE provided by
it. We consider the main neutrino-matter interactions (Just et al., 2015)
contributing to the neutrino energy deposition, $Q_{\mathrm{\nu}}$, which is
critical for triggering the explosion. It is well known that one-dimensional
(1d) spherically symmetric simulations do not explode due to the lack of
convection and hydrodynamical instabilities (e.g., Janka, 2012). Therefore, we
add an extra heating term in the gain region to launch explosions in 1d (see
e.g., Witt, 2020, for a similar approach).

We vary the setup only for $T\leq T_{\mathrm{th}}$, according to the following
different sets of models (see Table 3):

  1. i

One-dimensional simulations with no $\dot{E}_{\text{nuc}}$ and no neutrino
energy source term, $Q^{\mathrm{LD}}_{\mathrm{\nu}}$, in the regime
$T<T_{\mathrm{th}}$ that mainly corresponds to the infalling progenitor
shells. Note that we include neutrino heating at $T>T_{\mathrm{th}}$, where
the gain region is located, necessary to trigger the explosion. Four different
treatments of the composition have been employed: flashing scheme, RN16,
RN94eq148, and NSE. For the latter, we used the SFHo EOS also for
$T<T_{\mathrm{th}}$. In this way we check the impact of the composition in the
EOS, which should be negligible since it is dominated by the contributions of
electrons and photons.

  2. ii

One-dimensional simulations with $Q^{\mathrm{LD}}_{\mathrm{\nu}}$ and no
$\dot{E}_{\text{nuc}}$ for flashing scheme, RN16, RN94eq148, and NSE. This set
of models allows us to investigate the impact of neutrinos at low densities
within the different treatments and to examine the changes on the dynamics of
the explosions due to variations of the composition in the progenitor.

  3. iii

In order to show the impact of the nuclear energy generation, we perform one-
dimensional simulations adding $Q^{\mathrm{LD}}_{\mathrm{\nu}}$ and
$\dot{E}_{\text{nuc}}$ for the two different reduced networks. We consider
$\dot{E}_{\text{nuc}}$ in the two ways introduced in Section 2.2.2: everywhere
with $T<T_{\mathrm{inter}}$ and only below the shock. This makes it possible
to distinguish between the impact in the infalling layers and in the post-
shock region.

  4. iv

We run a total of three two-dimensional models, comparing the flashing scheme,
RN16, and RN94. For the three models we consider
$Q^{\mathrm{LD}}_{\mathrm{\nu}}$. For the latter two, $\dot{E}_{\text{nuc}}$
is taken into account.

We perform the 1d simulations using a grid with $n_{\mathrm{r}}=408$ zones
that are logarithmically spaced in radial direction with a central grid width
of $\Delta r=400\,\mathrm{m}$ and a maximum radius of $R_{\mathrm{out}}\simeq
9.05\times 10^{5}\,\mathrm{km}$. The 2d simulations run on a grid of
$n_{r}=400$ zones in the radial direction and $n_{\theta}=128$ in the angular
direction. We consider the employed resolution sufficient to adequately
resolve the hydrodynamic quantities (e.g., Obergaulinger & Aloy, 2017).

##  3 Results

###  3.1 Impact on the dynamics

In this section, we detail the impact that the different composition
treatments have on the dynamics of the explosion in the models introduced in
section 2.3.

Figure 2: Shock evolution of models without $\dot{E}_{\text{nuc}}$. Solid
lines stand for models with $Q^{\mathrm{LD}}_{\mathrm{\nu}}$, e.g. 1D_RN16.
Dotted lines correspond to models without, e.g. 1D_RN16*. The different colors
represent each treatment of the composition. The same colorcode is used in the
successive figures. Figure 3: Evolution of the ram pressure for the 1D_RN94
and 1D_RN16 models in the first $350\,\mathrm{ms}$ post-bounce. Bottom panel
shows the relative differences in the low density region for neutrino heating
(solid line) and shock radius (dotted). Figure 4: Mass shell evolution of
1D_RN16, 1D_RN16e and 1D_RN16E models, from left to right respectively. Grey
lines show the evolution of the mass shells from $1\,\mathrm{M_{\odot}}$ to
$2\,\mathrm{M_{\odot}}$ every $0.02\,\mathrm{M_{\odot}}$. The shock radius is
depicted in red. The inner green line corresponds to $T_{\mathrm{th}}$ and the
outer to $T_{\mathrm{inter}}$. Finally, $\dot{E}_{\mathrm{nuc}}$ is shown for
1D_RN16e and 1D_RN16E models. Figure 5: Shock evolution for the remaining
cases. Solid lines, again, correspond to the 1D_RN94 and 1D_RN16 models.
Dashed lines indicate for the 1D_RN16e and the 1D_RN94e models. Dash-dotted
lines correspond to the 1D_RN16E and the 1D_RN94E.

We start comparing 1D models with different composition treatment, no
$Q^{\mathrm{LD}}_{\mathrm{\nu}}$, and no $\dot{E}_{\mathrm{nuc}}$ at
$T<T_{\mathrm{th}}$: 1D_NSE*, 1D_flsh*, 1D_RN16*, and 1D_RN94*. Figure 2 shows
the evolution of the shock radius for these models (dotted lines). The
evolution for the latter three is identical, since their setup is effectively
very similar because of the absence of relevant energy source terms at low
temperatures. The main variation among these models is the composition, which
is an input for the Helmholtz EOS. We note that in the relevant temperature
regime the baryonic contribution is negligible and thus the way the nuclear
composition is treated has no dynamical effect. However, in the 1D_NSE* model,
small differences in the leptonic contribution to pressure and internal energy
can change the evolution of the shock.

On the other hand, in case $Q_{\mathrm{\nu}}$ is taken into account at
$T<T_{\mathrm{th}}$ (models 1D_flsh, 1D_NSE, 1D_RN16, and 1D_RN94, solid lines
in Figure 2), the evolution of the shock is now mainly influenced by this
additional energy source term, which is dominant over the changes between
different EOS. 1D_flsh and the 1D_NSE have similar behavior due to their
simplified composition, which leads to a similar amount of nucleons and, thus,
similar neutrino opacities and neutrino heating. Note that
$Q_{\mathrm{\nu}}\propto\kappa_{\mathrm{a}}\propto
n_{\mathrm{N}}\sigma_{\mathrm{\nu N}}$, where $\kappa_{\mathrm{a}}$ is the
absorption opacity, $n_{\mathrm{N}}$ is the free nucleon density, and
$\sigma_{\mathrm{\nu N}}$ is the neutrino-nucleon cross section.

Models 1D_RN16 and 1D_RN94 show a different evolution. The nuclei mapped into
simulations employing the RN94eq148 are different than the considered in the
original progenitor111 obtained in
https://2sn.org/nucleosynthesis/WH2007.shtml, which provides the ”aprox19”
composition. For this reason, RN94eq148 performs a readjustment where
$(\alpha,p)$ reactions are dominant in the oxygen rich layers. This leads to a
difference in the mass fractions of the neutrons of
$\Delta\mathrm{n}\equiv\log_{10}{(\frac{X^{\mathrm{RN94}}_{\mathrm{n}}}{X^{\mathrm{RN16}}_{\mathrm{n}}})}=0.3$
and $\Delta\mathrm{p}=2.3$ for protons, which which results an increase of the
opacity in the vicinity of the shock. This translates into differences of up
to $15\,\%$ in the neutrino heating in the low density regime (Figure 3)and in
a decrease of the ram pressure in the shock, which, although small, is
sufficient to modify its balance, allowing for an easier expansion in 1D_RN94.

Figure 6: Explosion energy of the models with no $\dot{E}_{\mathrm{nuc}}$
(solid), with $\dot{E}_{\mathrm{nuc}}$ below the shock (dashed), and with
$\dot{E}_{\mathrm{nuc}}$ in the low density region (dotted).

The energy from nuclear reactions significantly changes the dynamics. Figure 4
shows the evolution for simulations performed with the RN16 network taking
into account $Q^{\mathrm{LD}}_{\mathrm{\nu}}$ and the three different
configurations for $\dot{E}_{\mathrm{nuc}}$. The absence of
$\dot{E}_{\mathrm{nuc}}$ in the progenitor infalling shells in 1D_RN16 and
1D_RN16e lead to a very similar evolution of the shock and the mass shells.
The shock is revived after experiencing a fall back at $t\sim
300\,\mathrm{ms}$ and gaining enough energy from the system to expand.
Nevertheless, in 1D_RN16e, it expands slightly faster once the nuclear
reactions are starting to feed the system with energy (Figure 5). This is
better depicted in Figure 6, where the explosion energy shows that model
1D_RN16e is more than $15\,\mathrm{\%}$ more energetic than 1D_RN16 after
$500\,\mathrm{ms}$ post-bounce. On the other hand, model 1D_RN16E experiences
an early expansion with no fall-back and less layers are accreted to the PNS.
Initially, the Fe group rich layers absorb energy mainly through
$(\gamma,\alpha)$, $(n,\alpha)$, and $(\gamma,n)$ reactions, accelerating
slightly the collapse (e.g., Couch, 2017). The energy released in outer mass
shells, i.e. oxygen and silicon layers, is dominated by $(\alpha,\gamma)$ and
$(\mathrm{p},\gamma)$ reactions which produce $\dot{E}_{\text{nuc}}>0$. This
heats up the infalling matter substantially (see Figure 7). The contribution
of the nuclear reactions is larger than the one of neutrino absorption, for
$T\leq T_{\mathrm{th}}$, from $t=0.1\,\mathrm{s}$ on. Therefore, nuclear
reactions supply a significant amount of energy to the infalling layers, more
than $10^{50}\,\mathrm{erg\,s^{-1}}$, which is comparable to the change in
internal energy ($\sim 10^{50}-10^{51}\,\mathrm{erg\,s^{-1}}$) in the region.
This becomes the dominant heating source for $T\leq T_{\mathrm{th}}$,
increasing further the explosion energy. The additional energy source leads to
a decrease of the ram pressure in the shock and allows it to expand easier, in
agreement with Nakamura et al. (2014); Bruenn et al. (2006).

Figure 7: The upper panel shows the evolution of the integrated energy source
terms in the low-density region ($\mathcal{Q}$) for the models 1D_RN94E and
1D_RN16E in red and yellow, respectively.
$\mathcal{Q}^{\mathrm{LD}}_{\mathrm{nu}}=\int_{\mathrm{LD}}{Q^{\mathrm{LD}}_{\mathrm{nu}}\,dV}$
and
$\dot{\mathcal{E}}_{\mathrm{nuc}}=\int_{\mathrm{LD}}{\dot{E}_{\mathrm{nuc}}\,dV}$.
The lower panel shows the relative difference in the energy released by
nuclear reactions in the low density region. Figure 8:  Evolution of the
explosion energy (solid) and the nuclear energy generation (dotted) of the 2D
models (upper panel) and the shock radius (lower panel). Figure 9: Slices of
entropy and radial velocity in the 2D models at $500\,\mathrm{ms}$,
$1000\,\mathrm{ms}$ and $1500\,\mathrm{ms}$ of simulation. Figure 10:
Snapshots of the magnitude of $\dot{E}_{\text{nuc}}$ and $Y_{\mathrm{e}}$ at
of model 2D_RN94E.

Despite the similar overall evolution for 1D_RN16E and 1D_RN94E, we encounter
differences in the nuclear energy production. Before the bounce, as pointed
out previously, the photodisintegration reactions in the iron group absorb
energy from the environment. The RN16 (c.f., Table 1) is not able to track
efficiently these nuclei due to its simplicity. In contrast, the RN94eq148
(c.f., Table 2) has a broader set of nuclei in the relevant region of the
nuclear chart and, therefore, it is able to reproduce the released energy much
better. At intermediate times, between the bounce and $t\sim
500\,\mathrm{ms}$, the energy production in both networks is comparable as
shown by the fluctuations of $\Delta\dot{\mathcal{E}}_{\mathrm{nuc}}$ in
Figure 7. At these times, the $(\mathrm{p},\gamma)$ reactions behind the shock
are the main contributors to the energy generation. Finally, at latest times,
the differences stabilise around $50\,\mathrm{\%}$, where $(\alpha,\gamma)$
reactions become again dominant. A more detailed analysis of the relative
importance of the reactions is beyond the scope of this work and may be
explored in a further study.

All of the two-dimensional models show typical supernova explosions, with
explosion energies (Figure 8) of several $10^{50}\,\mathrm{erg}$ and, due to
axisymmetry, a prolate shape (Bruenn et al., 2016; Summa et al., 2016;
Vartanyan et al., 2018).

The nuclear energy generation modifies the explosion energy of the 2d models
in a consistent manner with the spherically symmetric models described above,
and leads to characteristic changes in the morphology (Figure 9). In 2D,
$\dot{E}_{\mathrm{nuc}}$ is less determinant for the onset of the explosion
than in 1D. The break of spherical symmetry allows for non-radial deformations
and convection which trigger an easier shock expansion. Actually, in 2D_RN16E,
the effect of $\dot{E}_{\mathrm{nuc}}$ on $R_{\mathrm{shock}}$ is small and
contributes to accelerate it around $t\sim 800\,\mathrm{ms}$ after bounce. At
this time, we notice that the evolution with respect to 2D_flsh is quite
different. There is an $\mathrm{\alpha}$-rich $\mathrm{\nu}$-driven wind
towards the southern pole (the cone filled with high-entropy, high-velocity
gas in Figure 9e). The rate at which nuclear reactions deposit energy in the
wind and post-shock region is significant, comparable to the growth rate of
$E_{\mathrm{exp}}$ (Figure 8). This supports the outflows and leads to a more
energetic explosion. It is interesting to see that this is not the case for
2D_RN94E. In this model, $\dot{E}_{\mathrm{nuc}}$ affects the evolution of the
shock at earlier times by decreasing the accretion of the infalling progenitor
material. In agreement with the 1D models, this helps the expansion and
favours less massive PNSs: $M_{\mathrm{PNS}}=1.67,\,1.72,$ and
$1.71\,\mathrm{M_{\odot}}$ for 2D_RN94E, 2D_flsh, and 2D_RN16E respectively.

At later times, downflows accreting onto the PNS squeeze its polar region,
launching high velocity, neutron-rich outflows ($Y_{\mathrm{e}}\sim 0.35$)
ejected from the vicinity of the PNS towards mainly the z direction. This is
shown in Figure 10, where we observe a neutron-rich bubble which moves from
the center (panel a) to $z\sim-6000\,\mathrm{km}$ in only $150\,\mathrm{ms}$
(panel b), and continues evolving towards higher latitudes mixing with the
surrounding material (panel c). The high $v_{\mathrm{r}}$, up to $\sim
0.16\,\mathrm{cm\,s^{-1}}$, avoid excessive neutrino absorption that would,
for the conditions prevailing here, increase the electron fraction. These
outflows tend to be amplified by the nuclear energy generation in the already
shocked areas. In model 2D_RN94E, their presence is not restricted to the
polar region only, but they are also emitted in the equatorial direction from
$t\sim 1\,\mathrm{s}$ of simulation (see Figure 10 as example). This geometry
may be due to the stochasticity involved in 2D simulations. The equatorial
outflow is sustained in time thanks to the nuclear energy released in the
region, $\dot{E}_{\mathrm{nuc}}\sim
10^{20}\,\mathrm{erg}\;\mathrm{g}^{-1}\mathrm{s}^{-1}$. As we show in the next
section, this mechanism allows for the production of more neutron rich
species.

###  3.2 Impact on the composition

Figure 11:  Integrated final ejecta composition of the 1D models. Red lines
correspond to the composition obtained in post-processing with the full
network WinNet in the 1D_RN16E (left) and 1D_RN94E (right). The post-
processing results for the 1D_flsh are displayed in grey for comparison. Green
dots stand for the values obtained from the network in situ, i.e. evolved in
the simulation. The values obtained with the same reduced network in post-
processing are depicted by orange diamonds. Figure 12: Integrated final
composition of the 2D models. Analogously from Figure 11, the grey line
corresponds to the post-processing calculation of model 2D_flsh.

We have studied the impact that reduced networks have on the composition of
the ejecta at the end of the simulation, i.e., $t=1.5\,\mathrm{s}$, in the
one- and two-dimensional models. First, we show the impact on abundances of
taking into account the nuclear energy generation when computing the post-
processing nucleosynthesis with WinNet.

In spherical symmetry, $\mathrm{Y_{e}}$ is very close to $0.5$. Therefore,
species beyond the iron group are not synthesized (Figure 11). This also
explains the only minor differences in the post-processing final composition
with WinNet for the 1D_flsh, the 1D_RN16E, and 1D_RN94E models, despite the
different dynamic evolution when the energy generation from the nuclear
reactions is taken into account.

In the two-dimensional models (Figure 12) a larger amount of nuclei are
involved due to the broader range of $\mathrm{Y_{e}}$. The 2D_flsh and
2D_RN16E model show standard CCSN nucleosynthesis, where mainly iron group
elements are formed together with lighter heavy species around $A\sim 90$ (see
e.g., Eichler et al., 2017; Harris et al., 2017; Wanajo et al., 2018; Witt et
al., 2021). In contrast, model 2D_RN94E produces larger amounts of heavier
species. The composition obtained with WinNet shows a significant enhancement
around the first r-process peak, $A\sim 80$, particularly
$\mathrm{{}^{84}Ge}$. Also, a small amount is observed up to the second
r-process peak, $A\sim 130$. This larger production of heavier species takes
place in the late low-$\mathrm{Y_{e}}$ outflows, which the energy released in
nuclear reactions helps to sustain. We observe an agreement between the post-
processing RN94eq148 (ex situ) and WinNet abundances in the 2D_RN94E model
until $A\sim 90$, showing that RN94eq148 is able to reproduce the main nuclei
synthesized in CCSN. There is a difference in $\mathrm{{}^{92}Mo}$ that can be
explained by the fact that it acts as a bottle neck for heavier species in
RN94eq148.

Figure 13: Chart with the RN94eq148 isotopes in boxes. Orange edges indicate
unstable nuclei and black stable ones. Bottom half of boxes depict in situ
integrated mass fractions for 2D_RN94E at the end of the simulation. Upper
half show the differences with respect to ex situ mass fraction, defined as
$\Delta X=\frac{X_{\mathrm{in}}}{X_{\mathrm{ex}}}$, for species with
$X_{\mathrm{i}}>10^{-5}$.

We show the composition obtained in situ in the models with the networks and
its respective post-processing, ex situ calculation. The goal is to analyse
which differences come from evolving the network together with the hydro. The
1D_RN16E and 1D_RN94E models show very good agreement in the most abundant
species, like $\mathrm{{}^{56}Ni}$. However, there are some discrepancies on
$\mathrm{{}^{30}S}$, $\mathrm{{}^{34}Ar}$, $\mathrm{{}^{40}Ca}$,
$\mathrm{{}^{44}Ti}$, $\mathrm{{}^{52}Cr}$, and $\mathrm{{}^{62}Zn}$ with
$\log_{10}(\Delta
X)\equiv\log_{10}{(\frac{X_{\mathrm{in}}}{X_{\mathrm{ex}}})}=0.5-1$. Lagragian
tracer particles, used to compute the composition in post-processing, are more
uncertain when tracking low density ejecta and especially the products of the
$\mathrm{\alpha}$-rich freeze out, as shown in Harris et al. (2017). This is
the origin of the systematic disagreement between in situ and ex situ
calculations in our models. In 1D, the representation of the trajectories is
much more simple, and accurate, reducing the uncertainty. The 2D_RN16E model
has differences of the same order as the previous ones in
$\mathrm{{}^{40}Ca}$, $\mathrm{{}^{44}Ti}$, $\mathrm{{}^{48}Cr}$, and
$\mathrm{{}^{52}Fe}$. The overproduction of such isotopes in situ is
consistent with the explosion morphology seen in Figure 9. It is characterized
by an extended high entropy and velocity shocked region where mainly
$\mathrm{\alpha}$-rich freeze out takes place. The in situ and ex situ
discrepancies in 2d_RN94E are depicted in Figure 13. It shows how the post-
processing underproduction of $\mathrm{{}^{44}Ti}$, $\mathrm{{}^{46}Ti}$,
$\mathrm{{}^{48}Ti}$, $\mathrm{{}^{48}Cr}$, $\mathrm{{}^{50}Cr}$,
$\mathrm{{}^{52}Cr}$, $\mathrm{{}^{54}Fe}$, and $\mathrm{{}^{56}Fe}$ lead to a
more neutron rich path with bottle necks in $\mathrm{{}^{80}Ge}$ and
$\mathrm{{}^{84}Se}$. In situ abundances show a more neutron-deficient pathway
starting from the aforementioned isotopes of $\mathrm{Ti}$, $\mathrm{Cr}$, and
$\mathrm{Fe}$, which produce more $\mathrm{{}^{92}Mo}$, $\log_{10}(\Delta
X_{\mathrm{{}^{92}Mo}})=1.35$.

##  4 Summary & Conclusions

We have presented a detailed study of how the treatment of the composition
within CCSN simulations impacts the explosion dynamics and nucleosynthesis.

We performed 1d and 2d CCSN simulations using the neutrino-hydrodynamics code
Aenus-Alcar (Just et al., 2015; Obergaulinger & Aloy, 2020). So far, this code
included the nuclear reactions outside the NSE regime only via the simplified
flashing scheme (Rampp & Janka, 2002), which gives neither accurate
information on the composition nor accounts for the energy release by nuclear
reactions. We used the reduced network module ReNet (see Appendix A) to
replace the flashing scheme by a $16\,\mathrm{\alpha}$-chain (RN16, see Table
1) and a $94$ isotopes network (RN94eq148, see Table 2). The latter is able to
reproduce the main nucleosynthesis yields in standard CCSN explosions (e.g,
Eichler et al., 2017). In addition, thanks to the $148$ nuclei considered in
steady state approximation, RN94eq148 is the most extended network in the
nuclear chart ever employed 222To date of submission. It should also be noted
that in terms of the number of isotopes, the one of the Oak Ridge group (e.g.,
Bruenn et al., 2020) is larger in state-of-the-art hydrodynamic simulations.
Both in-situ networks return the composition of the gas and the rate at which
nuclear reactions generate or consume internal energy.

The different composition in the low-density region have an impact on the
amount of nucleons, which can change the neutrino heating in the vicinity of
the shock. This modifies the ram pressure outside it and, therefore, its
evolution.

We have demonstrated how the energy released in the nuclear reactions impacts
the dynamics of the explosion. The energy generation in the pre-shocked
collapsing matter decreases, again, the ram pressure outside the shock and
allows it to expand easier, in agreement with Bruenn et al. (2006) & Nakamura
et al. (2014). On the other hand, the nuclear energy released in the shocked
region has a significant contribution, up to $20\,\%$, to the total explosion
energy. Differences between RN16 and RN94eq148 are small regarding the nuclear
energy generation, where $(\alpha,\gamma)$ and $(\mathrm{p},\gamma)$ are the
main channels of production. While the models presented are not very
energetic, we explored more energetic explosions, $E_{\mathrm{exp}}\sim
1\,\mathrm{B}$, and obtained similar impact.

Finally, we obtained the detailed nuclear yields of the models by applying the
nuclear network WinNet with $6545$ isotopes in an ex-situ, post-processing
step to Lagrangian tracer particles tracking the fluid flow. We compared its
results among different models and to the in-situ networks. In 1d, the
differences are small since the $Y_{\mathrm{e}}$ involved are very similar
among the models and close to $0.5$. In 2D, the variation in abundances among
different models get larger (Figure 12). The energy released in the nuclear
reactions helps to sustain late neutron-rich outflows ejected from the
vicinity of the PNS. The model 2D_RN94E shows how this mechanism allows weak
r-process to take place. Moreover, we have compared the final composition
obtained in situ and ex situ making use of RN16 and RN94eq148. In agreement
with Harris et al. (2017), we find significant discrepancies mainly in
products of the $\mathrm{\alpha}$-rich freeze out, since Lagrangian tracer
particles involve larger uncertainties when tracking such regions. Also, we
demonstrate how these differences can lead to variations on the
nucleosynthesis path which alter the final yields. Thus, it is necessary to
employ in situ realistic networks in CCSN simulations to obtain more accurate
ejecta composition.

We thank J. Austin Harris, W. Raphael Hix, Finia Jost and Andre Sieverding for
useful discussions. This work was supported by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) - Project-ID
279384907 - SFB 1245, the European Research Council under grant
EUROPIUM-667912, and the State of Hessen within the Research Cluster ELEMENTS
(Project ID 500/10.006). The authors gratefully acknowledge the computing time
provided to them on the high-performance computer Lichtenberg at the NHR
Centers NHR4CES at TU Darmstadt (Project-ID 1527). This is funded by the
Federal Ministry of Education and Research, and the state governments
participating. MR and MO acknowledge support by the Spanish Ministry of
Science, Education and Universities (PGC2018-095984-B-I00,
PID2021-127495NB-I00) and the Valencian Community (PROMETEU/2019/071). MO
acknowledges support from the Spanish Ministry of Science via the Ramón y
Cajal programme (RYC2018-024938-I).

## References

  * Anderson et al. (1999) Anderson, E., Bai, Z., Bischof, C., et al. 1999, LAPACK Users’ Guide, 3rd edn. (Philadelphia, PA: Society for Industrial and Applied Mathematics)
  * Balasi et al. (2015) Balasi, K., Langanke, K., & Martínez-Pinedo, G. 2015, Progress in Particle and Nuclear Physics, 85, 33, doi: 10.1016/j.ppnp.2015.08.001
  * Bethe (1990) Bethe, H. A. 1990, Rev. Mod. Phys., 62, 801, doi: 10.1103/RevModPhys.62.801
  * Bliss et al. (2020) Bliss, J., Arcones, A., Montes, F., & Pereira, J. 2020, Phys. Rev. C, 101, 055807, doi: 10.1103/PhysRevC.101.055807
  * Bruenn et al. (2006) Bruenn, S. W., Dirk, C. J., Mezzacappa, A., et al. 2006, in Journal of Physics Conference Series, Vol. 46, Journal of Physics Conference Series, 393–402, doi: 10.1088/1742-6596/46/1/054
  * Bruenn et al. (2016) Bruenn, S. W., Lentz, E. J., Hix, W. R., et al. 2016, ApJ, 818, 123, doi: 10.3847/0004-637X/818/2/123
  * Bruenn et al. (2020) Bruenn, S. W., Blondin, J. M., Hix, W. R., et al. 2020, ApJS, 248, 11, doi: 10.3847/1538-4365/ab7aff
  * Bugli et al. (2021) Bugli, M., Guilet, J., & Obergaulinger, M. 2021, MNRAS, 507, 443, doi: 10.1093/mnras/stab2161
  * Burrows et al. (2020) Burrows, A., Radice, D., Vartanyan, D., et al. 2020, MNRAS, 491, 2715, doi: 10.1093/mnras/stz3223
  * Burrows & Vartanyan (2021) Burrows, A., & Vartanyan, D. 2021, Nature, 589, 29, doi: 10.1038/s41586-020-03059-w
  * Cernohorsky & van Weert (1992) Cernohorsky, J., & van Weert, C. G. 1992, ApJ, 398, 190, doi: 10.1086/171848
  * Clayton (1968) Clayton, D. 1968, Principles of stellar evolution and nucleosynthesis: with a new preface (University of Chicago Press)
  * Couch (2017) Couch, S. M. 2017, Philosophical Transactions of the Royal Society of London Series A, 375, 20160271, doi: 10.1098/rsta.2016.0271
  * Couch et al. (2015) Couch, S. M., Chatzopoulos, E., Arnett, W. D., & Timmes, F. X. 2015, ApJ, 808, L21, doi: 10.1088/2041-8205/808/1/L21
  * Cowan et al. (2021) Cowan, J. J., Sneden, C., Lawler, J. E., et al. 2021, , 93, 015002, doi: 10.1103/RevModPhys.93.015002
  * Curtis et al. (2019) Curtis, S., Ebinger, K., Fröhlich, C., et al. 2019, ApJ, 870, 2, doi: 10.3847/1538-4357/aae7d2
  * Eichler et al. (2019) Eichler, M., Sayar, W., Arcones, A., & Rauscher, T. 2019, ApJ, 879, 47, doi: 10.3847/1538-4357/ab24cf
  * Eichler et al. (2017) Eichler, M., Nakamura, K., Takiwaki, T., et al. 2017, Journal of Physics G: Nuclear and Particle Physics, 45, 014001, doi: 10.1088/1361-6471/aa8891
  * Fields & Couch (2021) Fields, C. E., & Couch, S. M. 2021, ApJ, 921, 28, doi: 10.3847/1538-4357/ac24fb
  * Harris et al. (2017) Harris, J. A., Hix, W. R., Chertkow, M. A., et al. 2017, ApJ, 843, 2, doi: 10.3847/1538-4357/aa76de
  * Hix & Meyer (2006) Hix, W. R., & Meyer, B. S. 2006, Nuclear Physics A, 777, 188, doi: 10.1016/j.nuclphysa.2004.10.009
  * Hix & Thielemann (1999) Hix, W. R., & Thielemann, F. K. 1999, Journal of Computational and Applied Mathematics, 109, 321.  https://arxiv.org/abs/astro-ph/9906478
  * Intel (2016) Intel. 2016, Intel MKL PARDISO, https://software.intel.com/en-us/node/470282
  * Janka (2012) Janka, H.-T. 2012, Annual Review of Nuclear and Particle Science, 62, 407, doi: 10.1146/annurev-nucl-102711-094901
  * Janka et al. (2016) Janka, H.-T., Melson, T., & Summa, A. 2016, Annual Review of Nuclear and Particle Science, 66, 341, doi: 10.1146/annurev-nucl-102115-044747
  * Just et al. (2015) Just, O., Obergaulinger, M., & Janka, H. T. 2015, MNRAS, 453, 3386, doi: 10.1093/mnras/stv1892
  * Korobkin et al. (2012) Korobkin, O., Rosswog, S., Arcones, A., & Winteler, C. 2012, MNRAS, 426, 1940, doi: 10.1111/j.1365-2966.2012.21859.x
  * Kotake et al. (2006) Kotake, K., Sato, K., & Takahashi, K. 2006, Reports on Progress in Physics, 69, 971, doi: 10.1088/0034-4885/69/4/r03
  * Kuroda et al. (2020) Kuroda, T., Arcones, A., Takiwaki, T., & Kotake, K. 2020, ApJ, 896, 102, doi: 10.3847/1538-4357/ab9308
  * Lentz et al. (2015) Lentz, E. J., Bruenn, S. W., Hix, W. R., et al. 2015, ApJ, 807, L31, doi: 10.1088/2041-8205/807/2/L31
  * Lippuner & Roberts (2017) Lippuner, J., & Roberts, L. F. 2017, ApJS, 233, 18, doi: 10.3847/1538-4365/aa94cb
  * Martin et al. (2015) Martin, D., Perego, A., Arcones, A., et al. 2015, ApJ, 813, 2, doi: 10.1088/0004-637X/813/1/2
  * Mezzacappa et al. (2020) Mezzacappa, A., Endeve, E., Messer, O. E. B., & Bruenn, S. W. 2020, Living Reviews in Computational Astrophysics, 6, 4, doi: 10.1007/s41115-020-00010-8
  * Mösta et al. (2015) Mösta, P., Ott, C. D., Radice, D., et al. 2015, Nature, 528, 376, doi: 10.1038/nature15755
  * Müller (2020) Müller, B. 2020, Living Reviews in Computational Astrophysics, 6, 3, doi: 10.1007/s41115-020-0008-5
  * Müller et al. (2017) Müller, B., Melson, T., Heger, A., & Janka, H.-T. 2017, MNRAS, 472, 491, doi: 10.1093/mnras/stx1962
  * Müller (1986) Müller, E. 1986, A&A, 162, 103
  * Nakamura et al. (2022) Nakamura, K., Takiwaki, T., & Kotake, K. 2022, MNRAS, 514, 3941, doi: 10.1093/mnras/stac1586
  * Nakamura et al. (2014) Nakamura, K., Takiwaki, T., Kotake, K., & Nishimura, N. 2014, ApJ, 782, 91, doi: 10.1088/0004-637X/782/2/91
  * Obergaulinger & Aloy (2017) Obergaulinger, M., & Aloy, M. Á. 2017, MNRAS, 469, L43, doi: 10.1093/mnrasl/slx046
  * Obergaulinger & Aloy (2020) —. 2020, MNRAS, 492, 4613, doi: 10.1093/mnras/staa096
  * Obergaulinger et al. (2018) Obergaulinger, M., Just, O., & Aloy, M. A. 2018, Journal of Physics G Nuclear Physics, 45, 084001, doi: 10.1088/1361-6471/aac982
  * O’Connor & Couch (2018) O’Connor, E. P., & Couch, S. M. 2018, ApJ, 865, 81, doi: 10.3847/1538-4357/aadcf7
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3, doi: 10.1088/0067-0049/192/1/3
  * Rampp & Janka (2002) Rampp, M., & Janka, H. T. 2002, A&A, 396, 361, doi: 10.1051/0004-6361:20021398
  * Rauscher et al. (2002) Rauscher, T., Heger, A., Hoffman, R. D., & Woosley, S. E. 2002, ApJ, 576, 323, doi: 10.1086/341728
  * Reichert et al. (2022) Reichert, M., Obergaulinger, M., Aloy, M.-A., et al. 2022, arXiv e-prints, arXiv:2206.11914.  https://arxiv.org/abs/2206.11914
  * Reichert et al. (2021) Reichert, M., Obergaulinger, M., Eichler, M., Aloy, M. Á., & Arcones, A. 2021, MNRAS, 501, 5733, doi: 10.1093/mnras/stab029
  * Ristic et al. (2022) Ristic, M., Holmbeck, E. M., Wollaeger, R., et al. 2022, arXiv e-prints, arXiv:2206.02273.  https://arxiv.org/abs/2206.02273
  * Roberts et al. (2016) Roberts, L. F., Ott, C. D., Haas, R., et al. 2016, ApJ, 831, 98, doi: 10.3847/0004-637X/831/1/98
  * Sandoval et al. (2021) Sandoval, M. A., Hix, W. R., Messer, O. E. B., Lentz, E. J., & Harris, J. A. 2021, ApJ, 921, 113, doi: 10.3847/1538-4357/ac1d49
  * Schneider et al. (2019) Schneider, A. S., Roberts, L. F., Ott, C. D., & O'Connor, E. 2019, Phys. Rev. C, 100, doi: 10.1103/physrevc.100.055802
  * Sieverding et al. (2020) Sieverding, A., Müller, B., & Qian, Y. Z. 2020, ApJ, 904, 163, doi: 10.3847/1538-4357/abc61b
  * Steiner et al. (2013) Steiner, A. W., Hempel, M., & Fischer, T. 2013, ApJ, 774, 17, doi: 10.1088/0004-637X/774/1/17
  * Summa et al. (2016) Summa, A., Hanke, F., Janka, H.-T., et al. 2016, ApJ, 825, 6, doi: 10.3847/0004-637x/825/1/6
  * Timmes (1999) Timmes, F. X. 1999, ApJS, 124, 241, doi: 10.1086/313257
  * Timmes et al. (2000) Timmes, F. X., Hoffman, R. D., & Woosley, S. E. 2000, ApJS, 129, 377, doi: 10.1086/313407
  * Varma et al. (2022) Varma, V., Mueller, B., & Schneider, F. R. N. 2022, arXiv e-prints, arXiv:2204.11009.  https://arxiv.org/abs/2204.11009
  * Vartanyan et al. (2018) Vartanyan, D., Burrows, A., Radice, D., Skinner, M. A., & Dolence, J. 2018, MNRAS, 477, 3091, doi: 10.1093/mnras/sty809
  * Vartanyan et al. (2021) Vartanyan, D., Coleman, M. S. B., & Burrows, A. 2021, arXiv e-prints, arXiv:2109.10920.  https://arxiv.org/abs/2109.10920
  * Wanajo et al. (2018) Wanajo, S., Müller, B., Janka, H.-T., & Heger, A. 2018, ApJ, 852, 40, doi: 10.3847/1538-4357/aa9d97
  * Weaver et al. (1978) Weaver, T. A., Zimmerman, G. B., & Woosley, S. E. 1978, ApJ, 225, 1021, doi: 10.1086/156569
  * Winteler (2012) Winteler, C. 2012, PhD thesis, The University of Basel, Basel, Switzerland
  * Winteler et al. (2012) Winteler, C., Käppeli, R., Perego, A., et al. 2012, ApJ, 750, L22, doi: 10.1088/2041-8205/750/1/L22
  * Witt et al. (2021) Witt, M., Psaltis, A., Yasin, H., et al. 2021, ApJ, 921, 19, doi: 10.3847/1538-4357/ac1a6d
  * Witt (2020) Witt, M. R. 2020, PhD thesis, Technische Universität Darmstadt, Darmstadt, doi: https://doi.org/10.25534/tuprints-00011567
  * Wongwathanarat et al. (2017) Wongwathanarat, A., Janka, H.-T., Müller, E., Pllumbi, E., & Wanajo, S. 2017, ApJ, 842, 13, doi: 10.3847/1538-4357/aa72de
  * Woosley & Heger (2007) Woosley, S. E., & Heger, A. 2007, Phys. Rep., 442, 269, doi: 10.1016/j.physrep.2007.02.009
  * Woosley & Weaver (1995) Woosley, S. E., & Weaver, T. A. 1995, ApJS, 101, 181, doi: 10.1086/192237
  * Yasin et al. (2020) Yasin, H., Schäfer, S., Arcones, A., & Schwenk, A. 2020, Phys. Rev. Lett., 124, 092701, doi: 10.1103/PhysRevLett.124.092701
  * Yoshida et al. (2021) Yoshida, T., Takiwaki, T., Kotake, K., et al. 2021, ApJ, 908, 44, doi: 10.3847/1538-4357/abd3a3

##  Appendix A Reduced network implementation

We developed a highly flexible and adaptable reduced nuclear reaction network
code, ReNet. This reaction network solves the nuclear reaction network
equations (e.g., Clayton, 1968; Winteler, 2012; Hix & Thielemann, 1999; Hix &
Meyer, 2006; Lippuner & Roberts, 2017; Cowan et al., 2021):

$\displaystyle\dot{Y}_{i}$
$\displaystyle=\sum_{j}\mathcal{N}_{j}^{i}\lambda_{j}Y_{j}$ (Decays,
Photodisintegration)
$\displaystyle+\sum_{j,k}\frac{\mathcal{N}^{i}_{j,k}}{1+\delta_{jk}}\rho
N_{\mathrm{A}}\langle\sigma v\rangle_{j,k}Y_{j}Y_{k}$ (2-body)
$\displaystyle+\sum_{j,k,l}\frac{\mathcal{N}^{i}_{j,k,l}}{1+\Delta_{jkl}}\rho^{2}N_{\mathrm{A}}^{2}\langle\sigma
v\rangle_{j,k,l}Y_{j}Y_{k}Y_{l}.$ $\displaystyle\text{(3-body)},$ (A1)

with the abundance $Y_{i}$ of nucleus $i$, the amount of destroyed or created
nuclei within the reaction $\mathcal{N}$, the density $\rho$, the decay
constant $\lambda$, and the cross section $\langle\sigma v\rangle$. The
factors $\delta_{jk}$ and
$\Delta_{jkl}=\delta_{jk}+\delta_{kl}+\delta_{jl}+2\delta_{jkl}$ ensure that
there is no double counting of identical nuclei involved in the equation. This
set of coupled differential equations tend to be stiff due to the large range
of magnitudes of the reaction rates (e.g., Hix & Thielemann, 1999). The system
is therefore solved with an implicit Euler algorithm (e.g., Lippuner &
Roberts, 2017) to ensure numerical stability.

Every considered nucleus in Eq. (A) introduces an additional differential
equation. However, not all nuclei are connected with all others and the
resulting Jacobian is therefore sparse (e.g., Hix & Thielemann, 1999; Lippuner
& Roberts, 2017). We tested if the computational overhead of a sparse solver
(PARDISO, Intel 2016) is less than that of a dense matrix solver (LAPACK,
Anderson et al. 1999). We find a scaling of $\sim O(N^{1.6})$ and $\sim
O(N^{2.8})$ for a sparse and dense solver, respectively. The overhead of the
sparse solver makes this approach faster in case that we consider more than
$\sim 300$ nuclei (c.f., to $N\gtrsim 100$ in Winteler 2012). For the small
networks considered, we therefore use the dense LAPACK solver.

To reduce the set of considered nuclei without neglecting important reactions,
we assume a steady state for some chosen nuclei:

$\frac{\mathrm{d}Y_{i}}{\mathrm{d}t}=0.$ (A2)

The strategy of using a steady-state flow through intermediate nuclei was
already applied in other reduced networks333accessed via
https://cococubed.com/code_pages/burn_helium.shtml such as, e.g., the aprox13
(Timmes, 1999), aprox19 (Weaver et al., 1978), and aprox21 (Paxton et al.,
2011). As an example of this approximation, consider the three nuclei 24Mg,
28Si, and 27Al. 24Mg and 28Si are connected by an $\alpha$-capture. Another
possibility to synthesize 28Si from 24Mg is given by the reaction chain
${}^{24}\mathrm{Mg}(\alpha,\mathrm{p})^{27}\mathrm{Al}(\mathrm{p},\gamma)^{28}\mathrm{Si}$.
The set of differential equations (not showing protons and alphas) is given
by:

$\displaystyle\frac{\mathrm{d}Y_{\mathrm{{}^{24}Mg}}}{\mathrm{d}t}=$
$\displaystyle+\langle\sigma
v\rangle_{{}^{28}\mathrm{Si}(\gamma,\alpha)^{24}\mathrm{Mg}}\,Y_{{}^{28}\mathrm{Si}}$
(A3) $\displaystyle-\rho N_{\mathrm{A}}\langle\sigma
v\rangle_{{}^{24}\mathrm{Mg}(\alpha,\gamma)^{28}\mathrm{Si}}\,Y_{{}^{24}\mathrm{Mg}}\,Y_{\alpha}$
$\displaystyle-\rho N_{\mathrm{A}}\langle\sigma
v\rangle_{{}^{24}\mathrm{Mg}(\alpha,\mathrm{p})^{27}\mathrm{Al}}\,Y_{{}^{24}\mathrm{Mg}}\,Y_{\alpha}$
$\displaystyle+\rho N_{\mathrm{A}}\langle\sigma
v\rangle_{{}^{27}\mathrm{Al}(\mathrm{p},\alpha)^{24}\mathrm{Mg}}\,Y_{{}^{27}\mathrm{Al}}\,Y_{\mathrm{p}}$
$\displaystyle\frac{\mathrm{d}Y_{\mathrm{{}^{27}Al}}}{\mathrm{d}t}=$
$\displaystyle+\langle\sigma
v\rangle_{{}^{28}\mathrm{Si}(\gamma,\mathrm{p})^{27}\mathrm{Al}}\,Y_{{}^{28}\mathrm{Si}}$
(A4) $\displaystyle+\rho N_{\mathrm{A}}\langle\sigma
v\rangle_{{}^{24}\mathrm{Mg}(\alpha,\mathrm{p})^{27}\mathrm{Al}}\,Y_{{}^{24}\mathrm{Mg}}\,Y_{\alpha}$
$\displaystyle-\rho N_{\mathrm{A}}\langle\sigma
v\rangle_{{}^{27}\mathrm{Al}(\mathrm{p},\alpha)^{24}\mathrm{Mg}}\,Y_{{}^{27}\mathrm{Al}}\,Y_{\mathrm{p}}$
$\displaystyle-\rho N_{\mathrm{A}}\langle\sigma
v\rangle_{{}^{27}\mathrm{Al}(\mathrm{p},\gamma)^{28}\mathrm{Si}}\,Y_{{}^{27}\mathrm{Al}}\,Y_{\mathrm{p}}$
$\displaystyle\frac{\mathrm{d}Y_{\mathrm{{}^{28}Si}}}{\mathrm{d}t}=$
$\displaystyle-\langle\sigma
v\rangle_{{}^{28}\mathrm{Si}(\gamma,\mathrm{p})^{27}\mathrm{Al}}\,Y_{{}^{28}\mathrm{Si}}$
(A5) $\displaystyle-\langle\sigma
v\rangle_{{}^{28}\mathrm{Si}(\gamma,\alpha)^{24}\mathrm{Mg}}\,Y_{{}^{28}\mathrm{Si}}$
$\displaystyle+\rho N_{\mathrm{A}}\langle\sigma
v\rangle_{{}^{24}\mathrm{Mg}(\alpha,\gamma)^{28}\mathrm{Si}}\,Y_{{}^{24}\mathrm{Mg}}\,Y_{\alpha}$
$\displaystyle+\rho N_{\mathrm{A}}\langle\sigma
v\rangle_{{}^{27}\mathrm{Al}(\mathrm{p},\gamma)^{28}\mathrm{Si}}\,Y_{{}^{27}\mathrm{Al}}\,Y_{\mathrm{p}}.$

Setting $\frac{\mathrm{d}Y_{\mathrm{{}^{27}Al}}}{\mathrm{d}t}=0$, we obtain

$\displaystyle Y_{{}^{27}\mathrm{Al}}=$ $\displaystyle\frac{\langle\sigma
v\rangle_{{}^{28}\mathrm{Si}(\gamma,\mathrm{p})^{27}\mathrm{Al}}\,Y_{{}^{28}\mathrm{Si}}+\rho
N_{\mathrm{A}}\langle\sigma
v\rangle_{{}^{24}\mathrm{Mg}(\alpha,\mathrm{p})^{27}\mathrm{Al}}\,Y_{{}^{24}\mathrm{Mg}}\,Y_{\alpha}}{\rho\,N_{\mathrm{A}}\left(\langle\sigma
v\rangle_{{}^{27}\mathrm{Al}(\mathrm{p},\alpha)^{24}\mathrm{Mg}}+\langle\sigma
v\rangle_{{}^{27}\mathrm{Al}(\mathrm{p},\gamma)^{28}\mathrm{Si}}\right)Y_{\mathrm{p}}},$
(A6)

which can be inserted into Eq. (A3) and (A5), eliminating the necessary time
evolution of $Y_{\mathrm{{}^{27}Al}}$. We applied this procedure to a series
of nuclei using the Sympy Python package to derive the equations and the
Jacobian analytically. The latter is necessary because proton or neutron
abundances can end up in the denominator when we include all reactions that go
through the steady-state nucleus. In this way, we are able to evolve a large
region of the nuclear chart without explicitly evolving the abundance of all
nuclei.

We created a network with a minimum amount of nuclei, but still covering the
relevant region in the nuclear chart. Specifically, we focused on reproducing
the following quantities:

  * •

The nuclear energy $\dot{E}$.

  * •

Abundances of protons, neutrons, and alphas as well as other key nuclei such
as 44Ti and 56Ni.

  * •

Quantitative agreement of the abundance pattern at the end of the simulation.

In addition, we tested the reliability of ReNet. Therefore, we generated a
network with $\sim 800$ species (in the following referred to as ”full
network”), including all relevant nuclei up to Ge and reactions, however, not
including weak reactions444Because the nuclear flow evolves fairly close to
stability for CC-SNe and the simulation times are usually less than the half
life of the synthesized nuclei, this is a reasonable approximation. However,
for longer simulation times or very neutron-rich material that evolves far
from stability this is a point to improve.. This generated network is compared
to the well established reaction network WinNet (Winteler et al., 2012), which
was already used in previous publications (e.g., Winteler et al., 2012;
Korobkin et al., 2012; Martin et al., 2015; Eichler et al., 2019; Bliss et
al., 2020; Reichert et al., 2021; Ristic et al., 2022). For this, we use three
different typical explosive trajectories taken from the 2D_flsh model (see
Table 3). One trajectory is slightly neutron-rich ($Y_{e,5.8\,\mathrm{GK}}\sim
0.48$), one is symmetric ($Y_{e,5.8\,\mathrm{GK}}\sim 0.50$), and one is
slightly proton-rich ($Y_{e,5.8\,\mathrm{GK}}\sim 0.51$).

For WinNet and ReNet both, the nuclear energy as well as the final abundances
pattern are in excellent agreement (Fig. 14) as expected for an identical
input.

Figure 14: Left figure: Energy generation for two different trajectories.
Shown is the generated nuclear energy for four different network
architectures. The upper panel shows the energy generation for a slightly
proton-rich trajectory ($Y_{e,5.8\,\mathrm{GK}}\sim 0.51$), the middle panel a
trajectory with symmetric conditions and the lower panel a slightly neutron-
rich one ($Y_{e,5.8\,\mathrm{GK}}\sim 0.48$). The time is relative to the
start of the network calculation (i.e., $T=5.8\,\mathrm{GK}$). Right figure:
Mass fractions versus mass numbers at the end of the simulation for three
different trajectories and four different networks. The top panel shows the
result for a slightly proton-rich trajectory ($Y_{e,5.8\,\mathrm{GK}}\sim
0.51$), the middle panel a trajectory with symmetric conditions, and the lower
panel the result for a slightly neutron-rich one ($Y_{e,5.8\,\mathrm{GK}}\sim
0.48$).

For a slightly proton rich trajectory (i.e., $Y_{e,5.8\,\mathrm{GK}}\sim
0.51)$ the energy generation of all reaction networks quantitatively agree.
There are three episodes were the energy generation becomes significantly
lower or even negative for a full network (upper panel of Fig 14). This is
caused by an increase in temperature and, as as a consequence, a strong
negative energy contribution of $(\gamma,\mathrm{p})$ reactions. Even though
approximately implemented in RN94eq148, the energy generation differs
significantly in these episodes. The larger amount of considered nuclei of
RN94eq148 in comparison to RN16 does not improve the accuracy of the energy
generation since the differences are caused by photodissociation of nuclei
that are not included in RN94eq148. The final mass fractions on the other
hand, benefit from the larger amount of included nuclei (upper panel of Fig.
14). The differences between RN16 and RN94eq148 for nuclei with $28<A<52$ are
a direct consequence of additionally included nuclei such as 30Si or 30S,
located on the proton- and neutron- rich side off the diagonal in the nuclear
chart (see Figure 13). For the given trajectory, these nuclei are most
abundant. Since they are not included in RN16, the nuclear flow is forced to
the included symmetric nuclei. This causes 44Ti, but also other nuclei to be
overestimated by more than three magnitudes from the result of the full
network calculation. Given the fraction of around $6\%$ of ejected matter with
$Y_{e,5.8\,\mathrm{GK}}\geq 0.51$ in the 2D_flsh model, this can have a major
effect on 44Ti abundances at the end of the simulation. On the other hand,
44Ti is in excellent agreement ($0.9\%$) between RN94eq148 and the full
network. In addition, the amount of 56Ni differs for both reduced networks
only slightly from the full network ($1.5\%$ and $1.6\%$ for RN94eq148 and
RN16, respectively). Hereby the difference is even less than the one between
ReNet and WinNet ($3.5\%$) and of numerical origin.

The energy generation of a trajectory with symmetric conditions
($Y_{e,5.8\,\mathrm{GK}}=0.5$) evolves similarly to the previously discussed
proton-rich case. Again, the episodes of negative nuclear energy generation in
the full network are caused by $(\gamma,\mathrm{p})$ reactions on nuclei that
are neither included in RN16 nor in RN94eq148. In contrast to the slightly
proton-rich trajectory, the flow evolves along symmetric nuclei and the RN16
is therefore a better approximation compared to the previous, proton-rich
case. The agreement to the full network of 44Ti is with $3\%$ even better for
RN16 compared to the agreement of $35\%$ of RN94eq148. Similarly 56Ni is with
$10\%$ deviation more precise than for RN94eq148 with $22\%$ deviation. On the
other hand, protons are underrepresented in the composition of RN16. The
amount of nucleons is important for the neutrino opacities.

Even though only $6\%$ of the matter of the 2D_flsh run is neutron-rich
($Y_{e,5.8\mathrm{GK}}\leq 0.49$) it is still desirable to reproduce the
energy generation and composition within the reduced networks as a local
energy generation may also be able to change the dynamics of the simulation.
For these conditions, RN94eq148 is able to reproduce all major features of the
nuclear energy generation while RN16 differs by around one order of magnitude
(lower panel of Fig. 14). Furthermore, RN16 is unable to follow the episodes
of negative energy generation, because the nuclear flow is dominantly located
in the iron region which is not covered by RN16. Furthermore, due to the
mostly symmetric nuclei that are included in RN16, a lower electron fraction
can only be achieved by an increase of the neutron abundance.

Figure 15: Evolution of neutrons, protons, alphas, 44Ti, and 56Ni for a
neutron-rich trajectory. Shown are the results for four different reaction
networks.

As a consequence, the NSE composition of neutrons and protons already differs
significantly from the one of the full network (Fig. 15). We note that this
problem can be partly avoided by introducing another very neutron rich nuclei
(e.g., 56Cr within the aprox21, Paxton et al. 2011). The insufficient coverage
of the nuclear chart of RN16 leads also to a larger deviation of 44Ti compared
to RN94eq148 (factor $8$ instead of factor $5$ larger) as well as a larger
deviation of 56Ni ($88\%$ versus $41\%$). In general, the final abundances are
better reproduced with RN94eq148. Especially the peak around the most abundant
nucleus 60Ni (Fig. 14). However, the reduced network synthesizes even heavier
nuclei than the full reaction network. These heavy nuclei are even present
when calculating the same trajectory, including the steady state nuclei in the
network. Therefore, this is not an effect of the steady state approximation,
but rather of the exclusion of uneven nuclei in RN94eq148. This exclusion of
uneven nuclei in the calculation seem to bypass bottlenecks that are otherwise
present and that prevent the flow from synthesizing heavier elements.

Summarized, RN16 and RN94eq148 are good approximations to a full reaction
network that includes all nuclei. For more extreme conditions towards neutron-
and proton- rich conditions the larger reaction network RN94eq148 provides a
more accurate composition compared to the simpler RN16. The energy generation
is with the exception of neutron-rich conditions well reproduced for both
networks.


# A Topological Nomenclature for 3D Shape Analysis in Connectomics

Abhimanyu Talwar Zudi Lin Donglai Wei∗ Yuesong Wu Bowen Zheng
Jinglin Zhao Won-Don Jang Xueying Wang Jeff Lichtman Hanspeter Pfister
Harvard University
<EMAIL_ADDRESS>
###### Abstract

One of the essential tasks in connectomics is the morphology analysis of
neurons and organelles like mitochondria to shed light on their biological
properties. However, these biological objects often have tangled parts or
complex branching patterns, which make it hard to abstract, categorize, and
manipulate their morphology. In this paper, we develop a novel topological
nomenclature system to name these objects like the appellation for chemical
compounds to promote neuroscience analysis based on their skeletal structures.
We first convert the volumetric representation into the topology-preserving
reduced graph to untangle the objects. Next, we develop nomenclature rules for
pyramidal neurons and mitochondria from the reduced graph and finally learn
the feature embedding for shape manipulation. In ablation studies, we
quantitatively show that graphs generated by our proposed method align with
the perception of experts. On 3D shape retrieval and decomposition tasks, we
qualitatively demonstrate that the encoded topological nomenclature features
achieve better results than state-of-the-art shape descriptors. To advance
neuroscience, we will release a 3D segmentation dataset of mitochondria and
pyramidal neurons reconstructed from a 100$\mu$m cube electron microscopy
volume with their reduced graph and topological nomenclature annotations. Code
is publicly available at https://github.com/donglaiw/ibexHelper.

Figure 1: Illustration of topological nomenclature for mitochondria and
pyramidal neurons in EM connectomics. Given the input 3D segmentation, our
proposed algorithm generates reduced graph representation to capture essential
morphological structure and nomenclature for concise scientific communication.

##  1 Introduction

Recent advancements in large-scale electron microscopy (EM) allows the
generation of petabytes of serial images of brain tissue at nanometer
resolution [9, 20]. Machine learning methods have made automated 3D
reconstruction possible for individual neurons [7] and intracellular
organelles such as mitochondria [2]. Intriguingly, 3D shapes of these objects
resolved at the nano-scale are far more complicated than the classic depiction
in the neuroscience textbooks. Thus, novel morphology analysis tools are
required to advance our understanding of the basic properties of neuronal
compartments (Fig. 1).

Nevertheless, there are three main challenges. First, the branches and loops
of a non-convex object can tangle together in the 3D meshes, which make it
difficult for an intuitive perception of the underlying topology. Second,
there lacks an intuitive and concise way to convey the shape information of
neurons and organelles in the neuroscience community. Third, the goal of
traditional descriptors of 3D meshes is to compare objects with similar
scales, which is not suitable for the application on neurons and organelles
that have a wide range of spatial dimensions.

Figure 2: Overview of our topological nomenclature framework (first row) and
nomenclature rules (second row). Best view in color.

To tackle these challenges, we propose a topological nomenclature system to
abstract, categorize, and manipulate the 3D meshes of neurons and organelles
(Fig. 1). We first skeletonize them into vertexes and edges to untangle the
objects. We further prune them into a concise reduced-graph while preserving
the topological properties. To systematically name those objects, we propose a
nomenclature system borrowing ideas from nomenclature for organic compounds.
The primary aim of nomenclature in chemistry is to ensure that every name
refers to a specific compound without ambiguity. The naming systems, including
InChI [6] and SMILES, display more structural details but is more cumbersome
for scientific communication. Therefore in this work, we follow the IUPAC rule
[16] to generate the graph name that is more human-readable. To apply the
nomenclature system to shape analysis, we use the deep learning model for
self-supervised learning on graphs [11]. In comparison with traditional shape
descriptors like heat kernel signature (HKS) [14], one key characteristic of
our approach is that our graph representations are more intuitive to
understand than the HKS and allow for a simple shape decomposition into
primitives. Sundar _et al._ [19] explore the idea of the skeleton-based 3D
shape matching. They introduce the concept of a topological signature vector -
a low dimensional representation of a graph that can be a measure of
similarity. One difference from our approach is that their scheme generates an
acyclic skeletal graph that does not capture cycles or multiple paths between
two vertices.

To summarize, we present three main contributions in this paper. First, we
propose a shape abstraction method that converts 3D meshes into 2D graphs
using a novel combination of skeletonization and graph reduction to improve
the morphological perception. Second, by using the nomenclature system, we not
only make it more interpretable by neuroscientists but also further compress
the information needed for graph reconstruction. Third, we implement an
unsupervised model to embed these graphs into vector space for 3D shape
retrieval and decomposition.

##  2 Related Work

Chemical and Biological Nomenclature. Organic molecules, which contain carbon
as the backbone, exhibit a variety of structures. Therefore a concise and
intuitive nomenclature system is crucial for informative scientific
communication. The primary aim of nomenclature in chemistry and biology is to
ensure that every name refers to a specific compound without ambiguity. The
secondary aim is that the name can (to some extent) reflects a substance’s
structure. There are two main streams of chemical nomenclature. The first
stream that follows the IUPAC rule [16] is relatively simple and more human-
readable. Other systems, including InChI [6] and SMILES, display more
structural details in the name but is more cumbersome for scientific
communication. For large organic molecules like proteins that form a more
complex spatial arrangement, Flower [5] proposes to abstract the modules into
graphs and name them based on topology. Our work extends the nomenclature to
neurons and mitochondria morphologies in EM connectomics that do not have well
defined functional groups (e.g., benzene) and shape primitives (e.g.,
$\alpha-$ and $\beta-$helix of proteins) as biochemical compounds.

Figure 3: Illustration of graph reduction. We modify the breadth-first search
(BFS) traversal to consider only key nodes in the graph. Unlike BFS which
results in a tree, our traversal preserves loops within the graph to maintain
the cyclic topology of some mitochondria and pyramidal neurons that form
loops. Please see more details in Algorithm 1.

3D Shape Analysis in Connectomics. Recent advances in EM imaging have enabled
connectomics research at nano-scale resolution. For example, the 3D instance
masks of thousands of neurons and millions of intracellular organelles like
mitochondria are available for analysis. However, earlier studies only conduct
basic statistical analysis, including the length, volume, and surface-to-
volume ratio [9]. Recent Kanari _et al._ [8] propose a topological morphology
descriptor of neuronal structures based on distance transform to analyze the
shape of pyramidal neurons. However, such a method is not intuitive for human
perception as the reduced distance maps can hardly reflect the topology of the
original instances.

3D Shape Descriptors. Matching and retrieval of 3D shapes are mature
disciplines, and various successful schemes are out there. For example,
Ovsjanikov _et al._ [14] to create an isometry invariant shape retrieval
system by adopting a heat kernel signature (HKS) based deformation invariant
shape descriptor [18]. Sundar _et al._ [19] develop the skeleton-based 3D
shape matching algorithm and define a topological signature vector, which is a
low dimensional representation of a graph. These two methods are different
from our approach in two ways. First, our graph representations are more
understandable and intuitive compared to HKS while allowing a simple shape
decomposition into primitives. Unlike the existing skeleton-based 3D shape
matching scheme that generates an acyclic skeletal graph, our algorithm is
capable of capturing cycles or multiple paths between two vertices.

##  3 Method

In this section, we give a formal definition of our shape abstraction and
nomenclature system (Fig. 2). Given an input 3D mesh, we first transform it
into a reduced graph that preserves its topological information with a novel
skeletonization algorithm. We then determine its name in the nomenclature
system based on its category (e.g., mitochondria) and graph structure. We also
show how to compute the object feature in the nomenclature embedding space for
the following manipulation.

###  3.1 Topology-Aware Reduced Graph Generation

Starting from the volumetric representation of a 3D mesh, we convert it into a
reduced graph which preserves its topological structure, like the molecular
graph [12], using the following three steps:

Graph Initialization. We use an off-the-shelf skeletonization algorithm
proposed in Kálmán [15] to extract a 3D skeleton from the voxel
representation. We can view the extracted 3D skeleton as a weighted undirected
graph $G=(V,E,W)$ where $V\subset\mathbb{Z}^{*3}$ is the set of coordinates of
skeleton nodes in the 3D voxel grid, $E\subset V\times V$ is the set of edges,
and $W\subset\mathbb{R}^{+}$ is the set of edge weights.

Graph Reduction. Based on the degree of incident edges, we can divide the
skeleton vertices $V$ into junctions $J=\\{n\in V:degree(n)>2\\}$ and
endpoints $E=\\{n\in V:degree(n)=1\\}$. We aim to reduce the skeleton graph
$G$ to a graph $G_{S}$ whose set of vertices is $J\cup E$ (referred to as the
key nodes), and which preserves topological features of $G$ such as paths and
distances (along with the 3D skeleton) between any pair of key nodes. Further,
we also require $G_{S}$ to preserve any cycles present in the skeleton graph
$G$ and to preserve multiple paths between any two key nodes.

Figure 4: Illustration of extracting the longest path for the nomenclature assignment of undirected acyclic graphs (UAG) (Sec. 3.2). Although the longest-path problem in an arbitrary graph is $NP$-complete, the problem can be efficiently solved by running the breadth-first search (BFS) algorithm twice for UAGs. The first pass starts from an arbitrary leaf point (b), while the second pass starts from the farthest point found in the first round (c). Table 1: Prefixes and suffixes in our nomenclature system. We use Greek numeral prefixes to indicate the number of nodes on the longest path in an acyclic graph or on the rings, and use suffixes to categorize different type of instances. Prefixes (number of key points) | Suffixes (type)
---|---
1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | Mitochondira | Pyramidal Cell
mono- |  di- |  tri- |  tetra- |  penta- |  hexa- |  hepta- |  octo- |  ennea- |  deca- | -ito | -idal

For graph reduction, we modify the Breadth-First-Search traversal algorithm,
as outlined in Algorithm 1. At each step of the traversal, we only enqueue key
nodes to our traversal queue. We initialize the queue with any key node, and
while visiting a key node $v$, we only enqueue (1) any key nodes which are
adjacent to $v$, and (2) any other key nodes which are connected to $v$ by a
path in $G$ comprising only of non-key nodes (referred to as a simple path).
We define the “thickness” of an edge as the average of distance transforms of
its two vertices, and the “thickness” of a path can be calculated as the mean
of the thickness of each edge on the path weighted by edge lengths. During the
traversal, we keep track of two metrics for every pair of key nodes connected
by a simple path: (1) sum of lengths of all edges in $G$ on that simple path,
and (2) mean thickness of that simple path. An illustration of the graph
reduction scheme is shown in Fig. 3.

Graph Post-procession. In the biological systems, larger structures (_e.g_.,
large synapses with more vesicles and higher post-synaptic densities) usually
contribute more to the overall functionality. Therefore to further simplify
the graph without the distraction caused by numerous small structures, we can
remove edges and cycles whose path length is small relative to the total
length of edges in $G_{S}$. Thus, we contract all edges with a length lower
than a threshold value of $\tau$. With bigger $\tau$, we obtain a coarser-
level representation of the graph. In the experiments (Sec. 5), we perform
ablation studies to demonstrate how the threshold $\tau$ affect the quality of
the reduced graphs in terms of the agreement with human perception.

###  3.2 Topological Nomenclature Rules

Our nomenclature system for EM connectomics is modified upon the IUPAC
nomenclature of organic chemistry, which is not only invariant to the
deformation and the graph indexing order but also easily convertible back to
the graph representation. We add suffixes _-ito_ and _-idal_ to mitochondria
and pyramidal neurons, respectively.

Acyclic Graph. An acyclic graph is a graph having no cycles. A reduce graph
can be entirely a tree structure or contain tree branches. The nomenclature
rule for a tree is first to count the longest chain of vertexes, and assign a
prefix based on the number of vertexes. For example, the longest chain with
$n=5$ vertexes has a prefix _penta_ (Table 1). Find the longest path in a
general graph has been shown to be a $NP$-complete problem [3], but find the
longest path in an undirected tree graph can be solved efficiently by running
the breadth-first search (BFS) algorithm twice (see Fig. 4 for detail).
Therefore such a rule makes sure that for acyclic graphs, not only computer
programs but also human users can efficiently drive the corresponding
topological nomenclature. Then every vertex on the longest path is assigned a
location number from $1$ to $n$. For a branch, we use the location index as
prefix and name the branch recursively based on the rules. For simplicity, we
combine the prefix of branches with the same structure and omit the
description for branch topology if a branch only contains one node. An example
of the naming of the acyclic structure is shown in Fig. 2.

Algorithm 1 Skeleton Graph Traversal

1:procedure ReduceGraph

2: $Q\leftarrow\text{Any }v\in\left(J\cup P\right)$

3: $V_{S},E_{S},W_{S}\leftarrow\\{\\},\\{\\},\\{\\}$

4: $G_{S}\leftarrow(V_{S},E_{S},W_{S})$

5: $visited[v]\leftarrow\text{False }\forall v\in V$

6: while not Q.isEmpty() do

7: $src\leftarrow Q.dequeue()$

8: $visited[src]\leftarrow\text{True}$

9: $adj\\_list\leftarrow\text{KEY\\_NEIGHBORS}\left(G,src\right)$

10: for $trg\in adj\\_list$ do

11: if not visited[trg] then

12: $weight\leftarrow\text{PATH\\_LEN}(src,trg)$

13: if $(src,trg)\notin E_{S}$ then

14: $G_{S}\leftarrow\text{ADD\\_EDGE}(G_{S},(src,trg),weight)$

15: else$\triangleright$ True if multiple paths exist between $src$ and $trg$

16: $mid\leftarrow\text{NEW\\_KEY\\_NODE()}$

17: $G_{S}\leftarrow\text{ADD\\_EDGE}(G_{S},(src,mid),weight/2)$

18: $G_{S}\leftarrow\text{ADD\\_EDGE}(G_{S},(mid,trg),weight/2)$

19: end if

20: if $trg\notin Q$ then

21: $Q\leftarrow Q.enqueue(trg)$

22: end if

23: end if

24: if src == trg then $\triangleright$ True if $G$ has a cycle.

25: $mid1\leftarrow\text{NEW\\_KEY\\_NODE()}$

26: $mid2\leftarrow\text{NEW\\_KEY\\_NODE()}$

27: $G_{S}\leftarrow\text{ADD\\_EDGE}(G_{S},(src,mid1),weight/3)$

28: $G_{S}\leftarrow\text{ADD\\_EDGE}(G_{S},(mid1,mid2),weight/3)$

29: $G_{S}\leftarrow\text{ADD\\_EDGE}(G_{S},(mid2,trg),weight/3)$

30: end if

31: end for

32: end while

33:end procedure

Cyclic graph. If the reduced graph has circles, then we assign a higher
priority to the circles and name the graph accordingly. For the graph
structure with one circle, the prefix is _cyclo_. We then name branches use
parentheses containing relative location on the ring together with the branch
description described before. For a bicyclic graph where two circles share at
least one vertexes, the root numeral prefix of the graph name depends on the
total number of vertexes in all rings together[4]. The prefix _bicyclo_ denote
the sharing of at least two vertexes, while _spiro_ denote the sharing of only
one vertex. In between the prefix and the suffix, a pair of brackets with
numerals denotes the number of vertexes between each of the bridgehead ones.
These numbers are arranged in descending order and are separated by periods.
For example, a graph with a 3-vertex circle and a 5-vertex circle share two
vertexes (one edge) will be named $bicyclo[3,1,0]hexito$ (Fig. 2). Such rules
can be easily extrapolated into graphs with more than three circles. However,
in practice, we notice rare cases in the mitochondria and pyramidal neurons
with more than two circles.

###  3.3 Topological Nomenclature Embedding

In this subsection, we describe how to extract features from the reduced
graphs to estimate the commonality/dissimilarity between them. To represent a
graph as a matrix, we construct an adjacency matrix, whose elements are
connectivities between nodes in a graph. We employ a variational graph
autoencoder (VGAE) [11], which is a neural network for unsupervised learning
on graphs based on a variational autoencoder [10], to extract features for
each adjacency matrix.

We first normalize the adjacency matrix using the symmetric normalization
scheme, as done in VGAE [11]. Then, we feed the normalized adjacency matrix
into VGAE consisting of two graph convolutions and one fully connected layer
to reconstruct that matrix at the output. The dimensions of the two graph
convolutions are 32 and 16, respectively. We train the network for 200 epochs
to minimize the difference between the input adjacency matrix and the
reconstructed one and the Kullback-Leibler (KL) divergence of the embedding.
We use Adam optimizer with a learning rate of 0.01. We exploit the embedding
output of the graph convolutions as the nomenclature embedding.

Figure 5: Samples in our JWR-Mito300 and JWR-Pyr30 datasets. (a) Unlike
textbook illustrations, mitochondria can have complicated 3D shapes. (b)
Pyramidal neurons exhibit great diversity in the spatial distribution of
dendrites and their branching patterns.

##  4 Dataset

By inspecting the shape of pyramidal cells and mitochondria at nanometer
resolution, we found hundreds of structures that are very different from
textbook illustrations, which usually display simple spherical or tubular
structures. To further advance shape analysis in EM connectomics, we build a
dataset of non-trivial objects to exhibit the complexity of neuronal
structures and test our topological nomenclature system. We will release the
dataset publicly.

Data Acquisition. We imaged a tissue block from Layer II/III in the primary
visual cortex of an adult rat at a resolution of $4\times 4\times
30\leavevmode\nobreak\ nm^{3}$ using a multi-beam scanning electron microscope
(EM). After stitching and aligning the 2D images using multi-CPU clusters, we
obtained a final 3D image stack of 100 $\mu m$ cube.

3D Object Segmentation. Annotating the instances manually from scratch is not
feasible. However, the accuracy of existing automatic segmentation algorithms
can not generate object masks that are qualified enough for downstream
morphological analysis. To have a good tradeoff between segmentation quality
and efficiency, we first adopt the 3D U-Net model [17] for initial automatic
neuron and mitochondria segmentation. We then use a manual annotation tool [1]
to proofread and modify the segmentation results.

JWR-Mito300. We reconstructed all the mitochondria found in the somata of 11
cells: one pyramidal neuron, six interneurons, and four glial cells. Out of
thousands of mitochondira, we selected 316 of them that have nontrivial
topological structures with a volumetric size larger than 0.2 $\mu m^{3}$
(Fig. 5a).

JWR-Pyr30. For this dataset, we randomly selected 30 pyramidal cells whose
cell bodies are located in the central volume with the presence of a
significant portion (if not full) of their basal dendrites. Individual
pyramidal neurons have one apical dendrite pointing to the pial surface and an
axon often extending in the opposite direction. Nevertheless, they all show
distinct distributions of oblique and basal dendrites (Fig. 5b).

##  5 Experiments

In this section, we first quantitatively evaluate our nomenclature extraction
results in terms of the agreement with human perception on our JWR-Mito300 and
JWR-Pyr30 datasets. We then show qualitative results for two applications of
the extracted nomenclature feature, including 3D shape retrieval and 3D shape
decomposition.

Table 2: Ablation studies on graph reduction parameters on JWR-Mito300 and
JWR-Pyr30 dataset. We compute the average cosine similarity of the graph
Laplacian eigenvalues and the accuracy (correct if cosine similarity is bigger
than 0.95) for different design choices.

Dataset | Metric | BFS+Junction | BFS+Junction+Loop
---|---|---|---
$\tau$=0 |  $\tau$=4 |  $\tau$=0 |  $\tau$=4
JWR-Mito300 | Cosine | 0.891 | 0.949 | 0.914 | 0.956
Acc. | 0.376 | 0.785 | 0.493 | 0.854
JWR-Pyr30 | Cosine | 0.778 | 0.965 | N/A | N/A
Acc. | 0 | 0.767 | N/A | N/A

###  5.1 Nomenclature Extraction

Experimental Setup. As there is no known previous work that has similar
experiments, we conducted ablation studies on different design choices of our
proposed method. To generate the reduced graph, we started from the base
method of the modified BFS, referred to as BFS+Junction, and compared two
design choices. One is the edge length threshold $\tau$ for edge contraction
and the other is to preserve loop structure or not. Note that in the JWR-Pyr30
dataset, pyramidal neurons usually do not form loops and we only compare
different $\tau$ values.

Figure 6: Shape retrieval using graph features. Given (a) an input graph
representing the original 3D mesh, we show top-2 retrieval results using (b)
our topological nomenclature approach, (c) Heat Kernel Signature (HKS), and
(d) spectral embedding. The structures retrieved by our method are
qualitatively more similar to the query sample.

Evaluation Metric. To create ground truth nomenclature extraction labels for
JWR-Mito30 and JWR-Pyr30 datasets, we asked neuroscientists to draw their
perceived planar graph representation when showing them with the original 3D
object meshes. To quantitatively evaluate our automatic nomenclature
extractions, we use Cosine, the cosine similarity of the graph Laplacian
eigenvalues between the prediction and the ground truth, as the metric.
Empirically, we found that human labeling results have around 0.95 cosine
similarity due to the inherent ambiguity for small branches. Thus, we define a
prediction accuracy metric, Acc., with 0.95 cosine similarity as a threshold
for correctness.

Quantitative Results. As shown in Table 2, the choice of threshold $\tau$ for
small edge contractions is crucial, as the prediction accuracy almost doubles
with $\tau=4.0$ compared to $\tau=0$. Although all design choices have similar
average cosine similarity, they have a different number of correct predictions
that are acceptable for downstream analysis. Especially for pyramidal cells
from the JWR-Pyr30 dataset, the post-processing parameter $\tau=4$ helps to
remove many small spine structures.

Preservation of loops seems to have a small impact on overall performance, but
it is still essential for preserving the topology of cyclic mitochondria even
if they do not appear often. With the best design choices, our automatic
nomenclature extraction method can capture the essential topology of 3D
complex shapes without producing disturbing artifacts at around 80% accuracy.

Those results indicate that the graph extraction algorithm in the nomenclature
system produces high-quality representations that are consistent with human
perception. Considering that the nomenclature rule is designed to ensure that
every name refers to a specific structure without ambiguity, those results
further support the informativeness and conciseness of the assigned name for
scientific communication.

###  5.2 Shape Retrieval

Experimental Setup. For a given query 3D shape, users may want to find similar
shapes from the entire dataset. To this end, we perform 3D shape retrieval
using our topological nomenclature. The goal of this experiment is to find two
topological shapes that are most similar to a given query 3D shape. We use the
JWR-Mito300 dataset and randomly sample a query 3D shape from that dataset. We
discover two nearest neighbors of the query.

To compare two 3D shapes, we first measure pairwise differences between their
nomenclature embeddings by computing $L2$ distances. Then, we determine a
similarity between them as an average of matching costs. We use Hungarian
matching.

Qualitative Results. Fig. 6 shows 3D shape retrieval results of our algorithm
compared to both HKS [18] and spectral embedding [13]. The results exemplify
that our algorithm is capable of discovering topologically similar 3D shapes.
In contrast, HKS finds 3D shapes that have visually similar meshes but
different actual neuronal or mitochondria structures. Since the spectral
embedding encodes the entire graph, it fails to find relevant shapes.

Figure 7: 3D shape decomposition with topological nomenclatures. For
mitochondrion (top row) and a pyramidal neurons (bottom row), we show (a) the
learned part dictionary and (b) greedy decomposition result for an input
example.

###  5.3 Shape Decomposition Results

Experimental Setup. Decomposing topological nomenclatures into sub-graphs
enables users to understand the structures of 3D shapes. To define sub-graphs,
we construct a dictionary of our nomenclature embedding features.
Specifically, we apply $k$-means clustering algorithm to embedding features of
junctions (defined in Section 3.1) to generate words in the dictionary. We set
$k$ as 50 and 100 for the pyramidal neurons and mitochondria, respectively.
Note that we only use the junctions since end nodes have no local structures.
In the inference phase of decomposition, we perform matching between junctions
in a query nomenclature and the words in the dictionary. We first find a
junction with the minimum distance, and then remove it and its neighbor nodes
from the query nomenclature. We iterate this process until there are no more
junctions.

Qualitative Results. Fig. 7a visualizes dictionaries learned on the JWR-
Mito300 and JWR-Pyr30 datasets. It is observable that the words in the
dictionaries vary. Fig. 7b shows decomposition results of our nomenclatures.
Our method precisely decomposes the nomenclatures into sub-graphs.

##  6 Conclusion

In this paper, we introduced the topological nomenclature protocol to extract,
name, and manipulate the morphology of biological objects in EM connectomics.
We demonstrated the effectiveness of our novel nomenclature system through
quantitative ablation studies. Moreover, our unsupervised nomenclature
embedding successfully performed retrieval and decomposition of 3D shapes. We
will make the two datasets containing 316 mitochondria with complex morphology
and 30 pyramidal neurons publicly available. For future work, we will apply
our nomenclature scheme to a large-scale dataset to understand the diversity
and similarity of biological structures.

## References

  * [1] Daniel R Berger, H Sebastian Seung, and Jeff W Lichtman.  Vast: efficient manual and semi-automatic labeling of large 3d image stacks.  Frontiers in neural circuits, 12, 2018.
  * [2] Hsueh-Chien Cheng and Amitabh Varshney.  Volume segmentation using convolutional neural networks with limited training data.  In ICIP, 2017.
  * [3] Thomas H Cormen, Charles E Leiserson, Ronald L Rivest, and Clifford Stein.  Introduction to algorithms.  MIT press, 2009.
  * [4] Henri A Favre and Warren H Powell.  Nomenclature of organic chemistry: IUPAC recommendations and preferred names 2013.  Royal Society of Chemistry, 2013.
  * [5] Darren R Flower.  A topological nomenclature for protein structure.  Protein engineering, 11(9):723–727, 1998.
  * [6] Stephen Heller, Alan McNaught, Stephen Stein, Dmitrii Tchekhovskoi, and Igor Pletnev.  Inchi-the worldwide chemical structure identifier standard.  J. of cheminformatics, 5(1):7, 2013.
  * [7] Michał Januszewski, Jörgen Kornfeld, Peter H Li, Art Pope, Tim Blakely, Larry Lindsey, Jeremy Maitin-Shepard, Mike Tyka, Winfried Denk, and Viren Jain.  High-precision automated reconstruction of neurons with flood-filling networks.  Nature methods, 15(8):605, 2018.
  * [8] Lida Kanari, Paweł Dłotko, Martina Scolamiero, Ran Levi, Julian Shillcock, Kathryn Hess, and Henry Markram.  A topological representation of branching neuronal morphologies.  Neuroinformatics, 16(1):3–13, 2018.
  * [9] Narayanan Kasthuri, Kenneth Jeffrey Hayworth, Daniel Raimund Berger, Richard Lee Schalek, José Angel Conchello, Seymour Knowles-Barley, Dongil Lee, Amelio Vázquez-Reina, Verena Kaynig, Thouis Raymond Jones, et al.  Saturated reconstruction of a volume of neocortex.  Cell, 162(3):648–661, 2015.
  * [10] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  In ICLR, 2013.
  * [11] Thomas N Kipf and Max Welling.  Variational graph auto-encoders.  In NIPS Workshop on Bayesian Deep Learning, 2016.
  * [12] Alan D McNaught and Alan D McNaught.  Compendium of chemical terminology, volume 1669.  Blackwell Science Oxford, 1997.
  * [13] Andrew Y Ng, Michael I Jordan, and Yair Weiss.  On spectral clustering: Analysis and an algorithm.  In Advances in neural information processing systems, pages 849–856, 2002.
  * [14] M. Ovsjanikov, A. M. Bronstein, M. M. Bronstein, and L. J. Guibas.  A computer vision approach to isometry invariant shape retrieval.  In IEEE ICCV Workshops, 2009.
  * [15] Kálmán Palágyi.  A sequential 3d curve-thinning algorithm based on isthmuses.  In Advances in Visual Computing. Springer International Publishing, 2014.
  * [16] R Panico, WH Powell, and Jean-Claude Richer.  A guide to IUPAC Nomenclature of Organic Compounds.  Blackwell Scientific Publications, Oxford, 1993.
  * [17] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In MICCAI, 2015.
  * [18] Jian Sun, Maks Ovsjanikov, and Leonidas Guibas.  A concise and provably informative multi-scale signature based on heat diffusion.  Comput. Graph. Forum, 28, 07 2009.
  * [19] H. Sundar, D. Silver, N. Gagvani, and S. Dickinson.  Skeleton based shape matching and retrieval.  In Proceedings of the Shape Modeling International. IEEE Computer Society, 2003.
  * [20] Zhihao Zheng, J Scott Lauritzen, Eric Perlman, Camenzind G Robinson, Matthew Nichols, Daniel Milkie, Omar Torrens, John Price, Corey B Fisher, Nadiya Sharifi, et al.  A complete electron microscopy volume of the brain of adult drosophila melanogaster.  Cell, 174(3), 2018.