Core-collapse Supernova Simulations with Reduced Nucleosynthesis Networks

We present core-collapse supernova simulations including nuclear reaction networks that 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 number of nucleons and thus the neutrino-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 postshocked matter produces up to 20% more energetic explosions. 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, 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. These differences stem from the intrinsic characteristics of evolving the composition in hydrodynamic simulations or calculating it with Lagrangian tracer particles. The mass fractions of some Ca, Ti, Cr, and Fe isotopes are consistently underproduced in postprocessing calculations, leading to different nucleosynthesis paths. Our results suggest that large in situ nuclear reaction networks are important for a realistic feedback of the energy generation, the neutrino heating, and a more accurate ejecta composition.


INTRODUCTION
Massive stars (M ≳ 8 M ⊙ ) usually end their life in core-collapse supernovae (CCSN).These explosive events that are triggered by the collapse of the iron core and release several solar masses of products that were 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 a matter of study for many decades.In order to investigate the mechanisms involved, radiation-(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 neutrino 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 in 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 some several hundred nuclei (e.g., Woosley & Weaver 1995;Rauscher et al. 2002).Thus, including nuclear reaction networks in multidimensional CCSN 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 postprocessing (e.g., Winteler et al. 2012;Lippuner & Roberts 2017;Reichert et al. 2023b) using Lagrangian tracer particles.
If the temperatures are high enough (T ≳ 5 − 6 GK), fusion reactions and photodisintegrations are in socalled 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 (neglecting weak reactions).In NSE, the abundances are direct functions of the density, temperature, and electron fraction (Y e ) conditions.Thus, instead of the costly evolution of the reaction network, they can be used at high temperatures and be included in high-density EOS tables.If the temperature is lower, NSE no longer holds, and other methods are necessary to estimate the composition and energy generation.In this case, other simplified treatments are used.The socalled flashing scheme (Rampp & Janka 2002) takes into account neutrons, protons, α-particles, and a characteristic nucleus that depends on the thermodynamic conditions, considering instantaneous burning.For example, when silicon burning occurs, the silicon mass fraction (X( 28 Si)) is immediately transformed into 56 Ni and set to X( 28 Si) = 0.While the composition assumed by the flashing scheme is only a rough approximation, this scheme implicitly accounts for the generation or consumption of nuclear energy by the instantaneous burn-ing between the nuclei that is assumed to happen at the threshold temperature.Thus, explicit source terms for nuclear energy generation are not required.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, α-chains are the most commonly used for this purpose.
Several groups have employed reduced networks in their works (e.g., Bruenn et al. 2016Bruenn et al. , 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 αchain in the hydrodynamics, nuclear burning in the oxygen layer deposited additional pressure in the vicinity of the shock and assisted its expansion.Nakamura et al. (2014) performed simulations with a simplified lightbulb neutrino treatment using a 13 isotope α-network to study how the energy released by nuclear reactions affects the dynamics of the explosion.They concluded 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 ν-transport.They used a 14 α-chain and a 150 isotope network in axisymmetric (2D) models and studied the uncertainties derived from postprocessing nucleosynthesis.Their results showed the limitations of using postprocessing 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 affect the dynamics of the explosion and the nucleosynthesis outcomes.In Section 2 we introduce the radiationhydrodynamics code we used for the simulations (Section 2.1), the reduced network module implementation and 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.
Prior to the implementation of the reduced network module, Aenus-Alcar employed a nuclear EOS for matter at high densities (ρ > ρ th ∼ 10 7−8 g cm −3 ) assuming NSE between protons, neutrons, alpha particles, and a representative heavy nucleus.At lower densities, (ρ ≤ ρ th ), a Helmholtz-type EOS was used (see, e.g., Timmes & Arnett 1999;Timmes & Swesty 2000), 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 is 28 Si or 56 Ni, depending on the thermodynamic conditions.In this work, we have added a nuclear reaction network module to describe the composition outside of the NSE regime more accurately.

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 integration scheme in ReNet is a first-order implicit Euler scheme (as described in e.g., Hix & Thielemann 1999;Lippuner & Roberts 2017).The convergence criterion of the necessary Newton-Raphson scheme is given by the mass conservation, Y i A i = X i = 1, with the abundance Y i and mass number A i .This is identical to what is used in Lippuner & Roberts (2017).In contrast to postprocessing reaction networks, the thermodynamic quantities in the next time step are unknown for in-situ reaction networks.Therefore, ReNet burns hydrostatically in each time step with the temperature and density of the current step.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 we consequently use solvers for dense matrices.Within ReNet we furthermore assume some species to be in steady state in order to consider them indirectly and to avoid using a differential equation for them.This allows us to save computational cost (for more information, see Appendix A and see, 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 where Ėnuc is the variation in the internal energy due to the nuclear reactions, ∆m i is the mass excess of each species, and Ẏi is the abundance time variation of the nucleus i.

ReNet in Aenus-Alcar
As mentioned above, previous versions of Aenus-Alcar distinguish between two different density ranges where two different EOS regimes are used, one of which assumes NSE.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 this an inconsistency, we use the same temperature criteria in the NSE-network and the EOS transitions, such that above a certain temperature threshold (T th ), T > T th ∼ 5 − 6 GK, NSE is assumed both for the composition and the EOS.Otherwise, the 0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 network is employed together with the subnuclear EOS.Although 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 is of no concern.
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 th > T > T inter , where T inter is defined as the lower-limit temperature for the interpolation between EOSs (for a schematic representation, see Figure 2).Furthermore, at T = T th , the initial values for the network need to be given.These seeds are provided by an NSE solver that considers the same nuclei as the nuclear reaction network.When the temperature in a grid cell drops below T 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 ≤ T 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 ≤ T inter (regions II and III in Figure 2).On the other hand, the nuclear 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}.
energy generation is taken into account only after the shock has passed the cell, i.e., in the postshocked gas (region III).For this purpose, we use a simple shockdetection algorithm.We additionally consider the case without feedback of the nuclear energy generation from the nuclear reaction network.

Networks
In this work, we implemented two reaction network configurations into Aenus-Alcar: a 16-α chain (RN16), and a 94-species network (RN94).RN16 is similar to the widely used 19 isotope approximation network from Weaver et al. (1978), which approximately reproduces the nuclear energy generation within CCSN simulations and therefore provides feedback to the total energy of the system.RN94 is able to track the main species synthesized in the CCSNe (see 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 protonrich trajectories along stability up to 92 Mo.In addition, 148 species are considered in steady state, reproducing a network consisting of 242 nuclei.Therefore, at the time of writing, we present the most complete network 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 postprocessing with the full nuclear reaction network WinNet (Reichert et al. 2023b) and the reduced networks RN16 and RN94.We consider the unbound matter, i.e., matter with positive total energy and positive radial velocity, at the final time of the simulation (t f = 1.5 s).The evolution of the ejecta is followed by Lagrangian tracer particles calculated backward in time, as described in Reichert et al. (2023a) (see also Sieverding et al. 2022 for an in-depth analysis of the uncertainties that can arise with this method).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 M ⊙ , and each cell contains at least four tracers.At the end, the nucleosynthesis of each tracer particle is weighted with its corresponding mass.

Models
We have performed simulations (t f = 1.5 s) using the solar metallicity 20 M ⊙ mass progenitor from Woosley & Heger (2007).We map its precollapse composition to those of our two networks.We consider the EOS transition at T th = 5.8 GK1 , given that NSE breaks down around this temperature, and T inter = 5.0 GK.Since this work is focused on the impact of the composition treatment at T < T th , all of the models have the same configuration at high temperatures (areas I and IV in Figure 2).We apply the SFHo EOS (Steiner et al. 2013) and assume the composition provided by it in NSE.We consider the main neutrino-matter interactions (Just et al. 2015;Just et al. 2018) contributing to the neutrino energy deposition, Q ν , which is critical for triggering the explosion.It is well known that one-dimensional (1D) spherically symmetric simulations do not explode due 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.p.s indicates that it is switched on only in the postshocked region.The fourth column indicates whether the neutrino interactions are taken into account at ρ < 1.2•10 8 g cm −3 .Finally, the last column states for the dimensionality of the simulations.
to the lack of convection and hydrodynamical instabilities (e.g., Janka 2012).Therefore, we add additional heating factor, HF = 2.8, in the gain region to launch explosions in all of the 1D (see, e.g., Witt 2020, for a similar approach): , where Q gain ν is the neutrino energy deposition in the gain layer.In 2D models, we set HF = 1 and do not add any additional neutrino heating.In order to understand the impact of nuclear reaction networks coupled to the hydrodynamics, we have varied the treatment of the composition and nuclear energy generation at T ≤ T th .Four approaches have been used to describe the composition at low temperatures and densities, where a network is necessary: • The crudest description is to keep using the highdensity (or nuclear) EOS table at low temperatures as well.This is commonly done in many supernova simulations that focus on the early explosion phase (e.g., O'Connor & Couch 2018).The composition in EOS tables typically consists of neutrons, protons, alpha particles, and a representative nucleus.This nucleus is obtained from assuming NSE or from single nucleus approximation (SNA; see Schneider et al. 2017, for a discussion of the impact of these two treatments).These models are marked "SFHo" and their composition corresponds to NSE.
• The flashing scheme was the improved composition treatment that has been used in the previous studies with ALCAR.In this set of models, the version of the flashing scheme employed does not include the nuclear energy implicitly for comparison purposes.Therefore, to distinguish it from the original version (Rampp & Janka 2002), we refer to it as flashing † , and the models are labeled flsh.
• We use the two reduced networks described in Section 2.2.3.The model names are RN16 and RN94.
The composition obtained with the four approaches is different and influences the final abundances, but also the dynamics.Variations in the abundances of nuclei heavier than alphas have a minor contribution to the pressure because it is dominated by radiation and electrons, and the ion contribution is small.However, the changes in the number of neutrons and protons can impact the energy deposited by neutrinos.Therefore, we also ran a few toy models without the neutrino absorption at low densities, ρ < 1.2 • 10 8 g cm −3 (models with no Q LD ν ), i.e, Q ν = 0 in regions II, III, and IV of Figure 2.
Finally, we also used various models to explore the impact of the energy generation by nuclear reactions when using the two reduced networks in 1D and 2D.We consider Ėnuc in the two ways introduced in Section 2.2.2: everywhere with T < T inter (regions II and III in Figure 2) where the composition is determined by the reduced networks (models labeled E) and only between T < T inter and the shock, i.e., only region III (models labeled p.s.).We employ these two different Ėnuc configurations in order to distinguish between the impact in the infalling layers and in the postshock region.In addition, we ran a model with the flashing scheme including its nuclear energy generation (1D flshE) for comparison.
We perform the 1D simulations using a grid with n r = 408 zones that are logarithmically spaced in the radial direction with a central grid width of ∆r = 400 m and a maximum radius of R out ≃ 9.05×10 5 km.The 2D simulations run on a grid of n r = 400 zones in the radial direction and n θ = 128 in the angular direction.We consider the employed resolution sufficient to adequately resolve the hydrodynamic quantities (e.g., Obergaulinger & Aloy 2017).
The models computed are listed in Table 3.The overview of the results is given in Figure 3 where we present the shock and explosion energy evolution for all models.In the following sections, we discuss the impact of the composition and nuclear reactions on the dynamics and final abundances.

Impact on the dynamics
In this section, we study the effects that reduced networks have on the dynamics of the explosion.We show how the different treatments of the composition lead to changes in the neutrino absorption that modify the dynamics (Section 3.1.1)and the impact of the energy released by nuclear reactions on the dynamics of the explosion in 1D (Section 3.1.2)and 2D (Section 3.1.3)3.1.1.Composition First, we detail the impact that the different composition treatments have on the dynamics of the explosion in the models introduced in Section 2.3.
We start comparing 1D models with different composition treatment, no Q LD ν , and no Ėnuc at T < T th : 1D SFHo noQ LD ν , 1D flsh noQ LD ν , 1D RN16 noQ LD ν , and 1D RN94 noQ LD ν .The upper left panel of Figure 3 shows the evolution of the shock radius for these models (dotted lines).Without energy transfer between neutrinos and matter at low densities and temperatures, the main variation in the non-SFHo 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 in which the nuclear composition is treated has no dynamical effect.Thus, the composition differences between the non-NSE models do not affect the evolution of the shock.However, 1D SFHo noQ LD ν uses a different table for the leptonic contribution than the other three runs.Its slightly different pressure and internal energy modify the evolution of the shock.
When Q ν is taken into account in the entire simulation domain (models 1D flsh, 1D SFHo, 1D RN16, and 1D RN94, solid lines in the upper left panel of Figure 3), the different evolution of the shock is now mainly influenced by this additional energy source term, which dominates the changes between different EOS.1D flsh and 1D SFHo have an almost identical behavior due to their very similar composition, which consists of nucleons, alphas (in case of 1D SFHo), and a representative nucleus.This leads to a comparable number of nucleons and, thus, similar neutrino opacities and neutrino heating.Note that   absorption opacity, n N is the free nucleon density, and σ νN is the neutrino-nucleon cross section.
Models 1D RN16 and 1D RN94 show changes in the evolution due to the different species included in the simulation, which changes the number of nucleons, and therefore, the neutrino absorption.The nuclei mapped into simulations employing the RN94 are different than those that are considered in the original progenitor2 , which provides the aprox19 composition.For this reason, RN94 performs a readjustment where (α, p) reactions are dominant in the oxygen-rich layers.This leads to a difference in the mass fractions of the neutrons of ∆n ≡ log 10 ( 3 and ∆p = 2.3 for protons3 , which results in an increase in 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 4) and into a decrease in the ram pressure in the shock, which, although small, is sufficient to modify its balance, allowing for an easier expansion in 1D RN94.

Nuclear energy generation in 1D
The energy from nuclear reactions significantly changes the dynamics.The middle left panel of Figure 3 shows the shock evolution for simulations performed with the RN16 and RN94 networks, taking into account the three different configurations for Ėnuc .The impact of Ėnuc on the dynamics depends on whether the energy is released in the progenitor shells or the postshock region.To show it in a clear manner, we first focus on the different evolution in the models with RN16.The absence of Ėnuc in the progenitor infalling shells in 1D RN16 and 1D RN16e leads to a very similar evolution of the shock and the mass shells.The shock is revived after experiencing a fallback at t ∼ 300 ms and gains enough energy from the system to expand.Nevertheless, in 1D RN16e, it expands slightly faster when the nuclear reactions are starting to feed the system with energy.This is better depicted in the middle right panel of Figure 3, where the explosion energy shows that model 1D RN16e is more than 15% more energetic than 1D RN16 after 500 ms postbounce due to the nuclear energy released in the explosive nucleosynthesis.1D models with nuclear energy generation also in the progenitor shells (1D RN16E, 1D RN94E, and 1D flshE) experience an early expansion without fallback, and fewer layers are accreted to the PNS.Initially, the Fe-group-rich layers absorb energy mainly through (γ, α), (n, α), and (γ, n) reactions, which slightly accelerates the collapse (e.g., Couch 2017).The energy released in outer mass shells, i.e. oxygen and silicon layers, is dominated by (α, γ) and (p, γ) reactions, which produce Ėnuc > 0. This heats up the infalling matter substantially (see Figure 5).The contribution of the nuclear reactions is higher than that of neutrino absorption for T ≤ T th , from t = 0.1 s on.Therefore, nuclear reactions supply a significant amount of energy to the infalling layers, more than 10 50 erg s −1 , which is comparable to the change in internal energy (∼ 10 50 − 10 51 erg s −1 ) in the region.This becomes the dominant heating source for T ≤ T th , and further increases the explosion energy.The additional energy source leads to a decrease in the ram pressure on the shock and allows it to expand easier (Figure 4), in agreement with Nakamura et al. (2014); Bruenn et al. (2006).
We note that in terms of the shock radius and explosion energy, the evolution of model 1D flshE proceeds similarly to that of models 1D RN16E and 1D RN94E for the first ∼ 300 ms.During this period, the nuclear energy generation of the flashing scheme with its approximate composition of 28 Si and 56 Ni is not too different from that in the network models.At later times, reactions that cannot be represented by these two nuclei become more important and, consequently, the explosion energy of the flashing model starts to deviate from that of the network models (see the middle right panel of Figure 3).This behavior seems to indicate that a flashing scheme with more nuclei might reproduce the evolution of the network more closely.
Despite the similar overall evolution for 1D RN16E and 1D RN94E, we encounter some differences in the nu- clear energy production.Before the bounce, as pointed out previously, the photodisintegration reactions in the iron group absorb energy from the environment.The RN16 is not able to track these nuclei efficiently due to its simplicity.Because of this, we observe a spike of Ėnuc > 0, where (α, γ) reactions in the oxygen-and silicon-rich shells are able to release more energy than photodisintegrations absorb in the iron-and nickel-rich layer.In contrast, the RN94 has a broader set of nuclei in the relevant region of the nuclear chart and therefore includes a more complete set of reactions, such as e.g., additional (n, γ) reactions.This leads to the observed differences.At intermediate times, between the bounce and t ∼ 500 ms, 24 Mg(α, γ), 27 Al(p, γ), 28 Si(α, γ), and 54 Fe(p, γ) are the main contributors to the energy generation in the accreted layers.Because all these reactions are included in RN16 and RN94, the energy production in both networks at these times is comparable, as is shown by the fluctuations of ∆ Ėnuc in Figure 5. Finally, from t ∼ 500 ms on, the nuclear energy is mostly produced in the postshocked region.We observe how ∆ Ėnuc stabilizes around 50%.Part of this difference is due to the slightly different shock evolution, which evolves slightly faster in 1D RN16E.The remaining discrepancy cannot be easily broken down to an individual reaction, but rather to an overall slightly different nucleosynthetic path.The different number of included nuclei in the calculation causes the differences.While most of the main reactions contributing to Ėnuc in this region are included in both networks, e.g, 16 O(α, γ), 28 Si(α, γ), 54 Fe(p, γ), or 55 Co(p, γ), additional available paths in RN94 allow for the inclusion of several secondary reactions, which provide additional energy to the system, e.g, 58 Ni(p, γ) and 59 Cu(p, γ).Furthermore, a difference in the number of nucleons and α can lead to an additional deviation in the nuclear energy.

Nuclear energy generation in 2D
All of the 2D models show typical supernova explosions, with explosion energies (bottom panels of Figure 3) of several 10 50 erg and the characteristic prolate shape due to axisymmetry (Bruenn et al. 2016;Summa et al. 2016;Vartanyan et al. 2018).The nuclear energy generation also leads to changes in the morphology of the explosion in the 2D models (Figure 6).However, Ėnuc is less determinant for the onset of the explosion than in 1D because relaxing spherical symmetry allows for nonradial deformations and convection, which enhance the neutrino heating in the gain layer and trigger an easier shock expansion.
The long-term evolution, in particular the growth of the explosion energy, can differ quite significantly even between otherwise similar models.The shock waves of models 2D flsh and 2D RN16E expand anisotropically with a moderate north-south asymmetry at early times (left column in Figure 6).The comparable weakness of the shocks in one of the hemispheres allows for important accretion streams toward the PNS.These lead to a significant increase in the neutrino emission (and therefore in E exp ) that trigger an α-rich ν-driven outflow toward the southern pole (the cone filled with high-entropy high-velocity gas in Figure 6e) in 2D RN16E and an axial outflow in 2D flsh.In 2D RN16E, the effect of Ėnuc on R shock is small and contributes to accelerating it around t ∼ 800 ms after bounce.The rate at which nuclear reactions deposit energy in the wind and postshock region is significant, comparable to the growth rate of E exp (lower right panel of Figure 3).This supports the outflows and leads to a more energetic explosion.
The shock in 2D RN94E, in contrast, initially expands more symmetrically toward higher radii (bottom right panel of Figure 3), which we attribute to the higher amount of energy generated by nuclear reactions around t = 100 ms.This reduces the accretion inflows toward the PNS as a comparison of the masses indicates: M PNS = 1.67, 1.72, and 1.71 M ⊙ for 2D RN94E, 2D flsh, and 2D RN16E, respectively.Therefore, the number of neutrinos emitted after t ∼ 500 ms is significantly lower, which explains the absence of ν−driven winds in 2D RN94E and the stagnation of E exp , in contrast to 2D RN16E and 2D flsh.A is ejected in the polar direction and is eventually mixed up with external layers, B cannot overcome the accretion in the equatorial region and stagnates.Nevertheless, the high nuclear energy released in that region helps to sustain this outflow in time.
At later times, downflows accreting onto the PNS squeeze its polar region, launching high velocity neutron-rich outflows (Y e ∼ 0.35) ejected from the vicinity of the PNS.This is shown for model 2D RN94E in Figure 7, where we observe a neutron-rich bubble (indicated by the circle A) that moves from the center (panel a) to z ∼ −6000 km in only 150 ms (panel b), and continues evolving toward higher latitudes, mixing with the surrounding material (panel c).The high v r , up to ∼ 0.16 c, avoid excessive neutrino absorption that would increase the electron fraction for the conditions prevailing here.These outflows tend to be amplified by the nuclear energy generation in the already shocked areas.The predominantly polar direction of this bubble may at least partly be due to the assumed axisymmetry rather than being a general feature of such neutron-rich clumps.In fact, the presence of another such bubble (circle B in Figure 7) shows that they can be formed at any latitude as a consequence of the random dynamics of convection and accretion streams in the vicinity of the PNS.The lower radial velocity of the equatorial outflow (∼ 0.07 c) and the accretion of the upper layer prevent it from evolving toward a larger radius.Heating by nuclear energy reactions at rates around Ėnuc ∼ 10 20 erg g −1 s −1 plays an important role in sustaining the clump.As we show in the next section, this mechanism allows for the production of more neutron-rich species.
To sum up, in Section 3.1 we have studied the impact of reduced networks on the dynamics of the explosion in both 1D and 2D models.The different nuclei included can modify the number of nucleons and hence the neutrino absorption at T < T th .This different heat-ing in the infalling progenitor shells has an impact on the ram pressure and can therefore change the shock evolution.Analogously, the nuclear energy released in that region allows for a decrease of the ram pressure on the shock and favors an easier expansion.In contrast, in the shocked region, nuclear reactions are able to increase the explosion energy substantially.We demonstrated this by comparing the explosion energy of the models 1D RN16 and 1D RN94 ( Ėnuc = 0) to the explosion energy of the models 1D RN16e and 1D RN94e ( Ėnuc only in the shocked region).The increase in E exp when Ėnuc starts to release in the shocked region is striking for both cases (see the solid vs dashed lines in Figure 3).The 2D RN94 model suggests that the important nuclear energy generation in this region helps to sustain late low-Y e outflows in the equatorial direction.Finally, RN16 and RN94 show small differences in Ėnuc .

Impact on the nucleosynthesis
In this section, we study the effects that reduced networks have on the composition of the ejecta at the end of the simulation, i.e., t = 1.5 s, in 1D and 2D models.

Postprocessing with WinNet
First, we show the impact on the abundances by computing them with the postprocessing nucleosynthesis network WinNet.We compare the abundances obtained in models that include nuclear energy generation (1D RN16E, 1D RN94E, 2D RN16E, and 2D RN94E) with those from models that do not take Ėnuc into account (1D flsh and 2D flsh) in 1D and 2D, respectively.
In spherically symmetric models, Y e is very close to 0.5.Therefore, species beyond the iron group are not synthesized (Figure 8).This also explains the only minor differences in the postprocessing 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 2D models (Figure 9), a larger number of nuclei are involved due to the broader range of 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 ∼ 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 com-position obtained with WinNet shows a significant enhancement around the first r-process peak, A ∼ 80, particularly 84 Ge.Moreover, a small amount is observed up to the second r-process peak, A ∼ 130.This larger production of heavier species takes place in the late low-Y e outflows, which are supported by the nuclear energy release as outlined in Section 3.1.We observe an agreement between the postprocessing RN94 (ex situ) and WinNet abundances in the 2D RN94E model until A ∼ 90, showing that RN94 is able to reproduce the main nuclei that are synthesized in CCSN.There is a difference in 92 Mo that can be explained by the fact that it acts as a bottleneck for heavier species in RN94.Proton number 2.5 0.0 2.5 Log 10 ( X) Log 10 (X in ) Neutron number

Neutron number
Figure 10.Chart with the RN94 isotopes in boxes.The orange edges indicate unstable nuclei, and black edges show stable nuclei.The bottom half of the boxes depict in situ integrated mass fractions for 2D RN94E at the end of the simulation.The upper half show the differences with respect to the ex situ mass fraction, defined as ∆X = X in Xex , for species with Xi > 10 −5 .

Ex situ vs in situ calculations
In this section, we show the differences of computing the nucleosynthesis in the simulation (in situ) and in postprocessing (ex situ).For the comparison, we calculate the latter with the same network employed in the hydro (RN16 or RN94, depending on the model).
In Figures 8 and 9, we also show the composition obtained in situ in the models with the networks and its respective postprocessing, ex situ calculation.The goal is to analyze 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, such as 56 Ni.However, there are some discrepancies for 30 S, 34 Ar, 40 Ca, 44 Ti, 52 Cr, and 62 Zn with log 10 (∆X) ≡ log 10 ( Xin Xex ) ≈ 0.5 − 1. Lagrangian and Eulerian methods can both suffer from numerical errors, although in different ways.Lagrangian tracer particles, used to compute the composition in postprocessing, are more uncertain when tracking lowdensity ejecta and especially the products of the α-rich freeze-out, as shown in Harris et al. (2017).Furthermore, accurately representing features varying on short length and time scales such as those arising from multidimensional potentially turbulent flows can be difficult and require an exceedingly large number of tracers.The numerical diffusivity that Eulerian methods require for stability can lead to artificial mixing of species across the grid.While its importance decreases with increas-ing grid resolution, much finer grids than feasible are necessary to suppress this effect.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 simpler and accurate, reducing the uncertainty to that arising mainly from the aforementioned numerical diffusion.We therefore tested its impact based on the 1D RN94E model by recalculating it with different resolutions of 500, 600, 700, and 1000 radial zones.While the ex situ abundances for all models were almost identical, with maximum deviations of the order of ∼ 40% for 42 Ti and 54 Ni, some elements (i.e., 30 S, 34 Ar, and 62 Zn) showed variations in the order of ∼ 60% − 75% for the in-situ abundances.All other elements also converged within the in-situ networks.The impact of numerical diffusion could be decreased in further studies by modifying the advection scheme of the composition (Plewa & Müller 1999).
The 2D RN16E model has differences on the same order as the previous ones in 40 Ca, 44 Ti, 48 Cr, and 52 Fe.The overproduction of these isotopes in situ is consistent with the explosion morphology seen in Figure 6.It is characterized by an extended high-entropy and highvelocity shocked region in which mainly α-rich freezeout takes place.The in situ and ex situ discrepancies in 2D RN94E are depicted in Figure 10.Some differences in the overall nucleosynthesis trend can be observed by looking at the most neutron-rich and deficient isotopes of Zn, Ge, Se, Kr, Sr, and Zr.The postprocessing underproduction of 44 Ti, 46 Ti, 48 Ti, 48 Cr, 50 Cr, 52 Cr, 54 Fe, and 56 Fe predominantly leads to a more neutron-rich path with bottlenecks in 80 Ge and 84 Se.In situ abundances show a tendency of being more neutron-deficient, starting from the aforementioned isotopes of Ti, Cr, and Fe. Figure 10 shows that matter accumulates at the proton-rich isotopes 72 Se, 74 Se, and 76 Se, suggesting that even the highly extended RN94 network still has some bottlenecks that could be better represented with an even more extended network.Finally, this more neutron-deficient path ends up producing more 92 Mo, log 10 (∆X92 Mo ) = 1.35.

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;Just et al. 2018;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 assumes that the gas consists only of nucleons and a representative nucleus, for which, depending on the temperature, we use 28 Si or 56 Ni.We used the reduced network module ReNet (see Appendix A) to replace the flashing scheme by a 16 α-chain (RN16) and a 94-isotope network (RN94).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, RN94 is the most extended network in the nuclear chart ever employed4 in state-of-the-art hydrodynamic simulations.Both insitu networks return the composition of the gas and the rate at which nuclear reactions generate or consume internal energy.
The different compositions in the low-density region have an impact on the number of nucleons, which can change the neutrino heating in the vicinity of the shock.This modifies the ram pressure outside of 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 preshocked collapsing matter again decreases the ram pressure outside the shock and allows it to expand easier, in agreement with Bruenn et al. (2006) & Nakamura et al. (2014).The nuclear energy released in the shocked region has a significant contribution, up to 20%, to the total explosion energy.The flashing scheme with 28 Si and 56 Ni is not able to reproduce the nuclear energy generation in this region, which leads 1D flshE to a lower explosion energy than 1D RN16E and 1D RN94E.Nevertheless, their overall evolution is similar.The differences between RN16 and RN94 are small regarding the nuclear energy generation, where (α, γ) and (p, γ) are the main production channels.While the models we presented are not very energetic, we explored more energetic explosions, E exp ∼ 1 B, and obtained a similar impact.
Finally, we obtained the detailed nuclear yields of the models by applying the nuclear network WinNet with 6545 isotopes in an exsitu postprocessing step to Lagrangian tracer particles tracking the fluid flow.We compared the postprocessing results among different models and to the in-situ networks.In 1D, the differences are small because the Y e involved are very similar among the models and are close to 0.5.In 2D, the variation in abundances among different models becomes stronger (Figure 9).The energy released in the nuclear reactions helps to sustain late neutron-rich outflows ejected from the vicinity of the PNS.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 RN94.In agreement with Harris et al. (2017), we find significant discrepancies mainly in products of the α-rich freeze-out because Lagrangian tracer particles involve larger uncertainties when tracking such regions.For the in situ results, we identify the resolution-dependent numerical diffusion of species with low abundances as a factor contributing to the discrepancies.Moreover, we demonstrate how these uncertainties propagate, leading to variations in the nucleosynthesis path that alter the final yields.
What are the advantages and disadvantages of evolving a network in the simulations?While Lagrangian tracer methods and the ex situ results based on them are well suited for dense regions (e.g., Price & Federrath 2010) and avoid the excess numerical diffusion that may beset grid-based Eulerian schemes and, consequently, in situ abundances, they lack mixing and are more uncertain in tracking low-density regions and their nucleosynthesis, e.g. products of the α−rich freeze-out.This work suggests that it is necessary to employ in situ realistic networks in CCSN simulations with a fine grid resolution, or with a less numerically diffusive advection scheme, to obtain a realistic feedback of the energy generation, the neutrino opacities, and a more accurate ejecta composition.Thus, this study showed the strengths and weaknesses of employing networks in CCSN simulations and, hopefully, can help future simulations to decide depending on the goal of the study.in RN94, the energy generation differs significantly in these episodes.The larger amount of considered nuclei of RN94 in comparison to RN16 does not improve the accuracy of the energy generation because the differences are caused by photodissociation of nuclei that are not included in RN94.The final mass fractions benefit from the larger number of included nuclei (upper panel of Fig. 11).The differences between RN16 and RN94 for nuclei with 28 < A < 52 are a direct consequence of additionally included nuclei such as 30 Si or 30 S, located on the proton-and neutron-rich side of the diagonal in the nuclear chart (see Figure 10).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 44 Ti, 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 GK ≥ 0.51 in the 2D flsh model, this can have a major effect on 44 Ti abundances at the end of the simulation.On the other hand, 44 Ti is in excellent agreement (0.9%) between RN94 and the full network.In addition, the amount of 56 Ni differs only slightly for both reduced networks from the full network (1.5% and 1.6% for RN94 and RN16, respectively).Hereby, the difference is even smaller than the one between ReNet and WinNet (3.5%) and is of numerical origin.
The energy generation of a trajectory with symmetric conditions (Y e,5.8 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 (γ, p) reactions on nuclei that are neither included in RN16 nor in RN94.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 with the full network of 44 Ti is with 3% even better for RN16 than the agreement of 35% of RN94.Similarly, 56 Ni is with 10% deviation more precise than for RN94 with 22% deviation.On the other hand, protons are underrepresented in the composition of RN16.The number 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.8GK ≤ 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, RN94 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. 11).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.As a consequence, the NSE composition of neutrons and protons already differs significantly from that of the full network (Fig. 12).We note that this problem can be partly avoided by introducing another very neutron-rich nucleus (e.g., 56 Cr within the aprox21, Paxton et al. 2011).The insufficient coverage of the nuclear chart of RN16 also leads to a larger deviation of 44 Ti compared to RN94 (larger by a factor 8 instead of a factor 5) as well as a larger deviation of 56 Ni (88% versus 41%).In general, the final abundances are better reproduced with RN94, especially the peak around the most abundant nucleus 60 Ni (Fig. 11).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 RN94.This exclusion of uneven nuclei in the calculation seems to bypass bottlenecks that are otherwise present and that prevent the flow from synthesizing heavier elements.
To summarize, RN16 and RN94 are good approximations to a full reaction network that includes all nuclei.For more extreme conditions toward neutron-and proton-rich conditions the larger reaction network RN94 provides a more accurate composition than the simpler RN16.The energy generation is well reproduced for both networks, with the exception of the neutron-rich conditions.

Figure 1 .
Figure 1.Density and temperature achieved in a characteristic CCSN simulation.It focuses on the region of the ρ-T plane close to the transition between the two EOS.The color depicts the amount of mass that fulfilled each (ρ, 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.

Figure 2 .
Figure 2. Schematic representation of the different regions of interest in this work.It shows the mass shell evolution (gray) of a characteristic 1D model.The red line shows the shock evolution.The temperatures involved in the transition between the NSE and the network regimes are depicted in green.The purple line shows the evolution of ρ LDν , the density at which we switch off the neutrino absorption in the toy models described in Section 2.3.

Figure 3 .
Figure 3. Shock evolution (left panels) and explosion energy (right panels) of the models discussed in this work, summarized in Table3.The different colors represent each treatment of the composition.The solid lines in the upper and middle panels stand for 1D models without Ėnuc, e.g.1D RN16.The dotted lines correspond to models without Ėnuc or Q LD ν , e.g.1D RN16 noQ LD ν .The dashed lines indicate the 1D RN16e and the 1D RN94e models.The dash-dotted lines correspond to the 1D RN16E, the 1D RN94E, and the 1D flshE.The solid lines in the lower panels show the average shock radius (left) and explosion energy (right) for the 2D models.For the latter, we show the integrated nuclear energy generation Ėnuc = T<T th Ėnuc dV (dotted).

Figure 4 .
Figure 4. Evolution of the ram pressure for the 1D RN94, 1D RN16, and 1D RN16E models in the first 350 ms postbounce.The bottom panel shows the relative differences in the integrated neutrino energy deposition at T < T th , Q T<T th ν = T<T th Qν dV , (solid line) and shock radius (dotted) of 1D RN94 and 1D RN16.

Figure 5 .
Figure 5.The upper panel shows the evolution of the integrated energy source terms in the T < T th region (Q) for the models 1D RN94E and 1D RN16E in red and yellow, respectively.Q T<T th ν = T<T th Qν dV and Ėnuc = T<T th Ėnuc dV .The lower panel shows the relative difference in the energy released by nuclear reactions in the low-density region.

Figure 6 .
Figure6.Slices of entropy and radial velocity in the 2D models at 500 ms, 1000 ms, and 1500 ms of simulation.

Figure 7 .
Figure7.Snapshots of the magnitude of Ėnuc, and Ye of model 2D RN94E at t = 1100 ms, t = 1250 ms, and t = 1400 ms, respectively.We highlight the different evolution of two low-Ye clumps (A and B) ejected from the vicinity of the PNS.While A is ejected in the polar direction and is eventually mixed up with external layers, B cannot overcome the accretion in the equatorial region and stagnates.Nevertheless, the high nuclear energy released in that region helps to sustain this outflow in time.

Figure 8 .Figure 9 .
Figure 8. Integrated final ejecta composition of the 1D models.The red lines correspond to the composition obtained in postprocessing with the full network WinNet in the 1D RN16E (left) and 1D RN94E (right).The postprocessing results for the 1D flsh are displayed in gray for comparison.The 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 postprocessing are depicted by orange diamonds.

Figure 11 .
Figure11.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 (Ye,5.8GK ∼ 0.51), the middle panel a trajectory with symmetric conditions and the lower panel a slightly neutron-rich one (Ye,5.8GK ∼ 0.48).The time is relative to the start of the network calculation (i.e., T = 5.8 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 (Ye,5.8GK ∼ 0.51), the middle panel a trajectory with symmetric conditions, and the lower panel the result for a slightly neutron-rich one (Ye,5.8GK ∼ 0.48).

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

Table 3 .
List of the different models