Implications for Post-Processing Nucleosynthesis of Core-Collapse Supernova Models with Lagrangian Particles

TL;DR

This paper examines the uncertainties in post-processing nucleosynthesis calculations in 2D core-collapse supernova models, highlighting limitations of tracer particles and proposing detailed analysis of these effects across multiple progenitor models.

Contribution

It provides a comprehensive analysis of uncertainties in nucleosynthesis post-processing in 2D CCSN models using tracer particles, across different progenitor masses.

Findings

Tracer particle resolution affects nucleosynthesis accuracy.

Thermodynamic inconsistencies impact composition evolution.

Uncertainties vary with progenitor mass and explosion dynamics.

Abstract

We investigate core-collapse supernova (CCSN) nucleosynthesis with self-consistent, axisymmetric (2D) simulations performed using the radiation-hydrodynamics code Chimera. Computational costs have traditionally constrained the evolution of the nuclear composition within multidimensional CCSN models to, at best, a 14-species $\alpha$-network capable of tracking only $(\alpha,\gamma)$ reactions from $^{4}$He to $^{60}$Zn. Such a simplified network limits the ability to accurately evolve detailed composition and neutronization or calculate the nuclear energy generation rate. Lagrangian tracer particles are commonly used to extend the nuclear network evolution by incorporating more realistic networks in post-processing nucleosynthesis calculations. However, limitations such as poor spatial resolution of the tracer particles, inconsistent thermodynamic evolution, including misestimation of expansion timescales, and uncertain determination of the multidimensional mass-cut at the end of the simulation impose uncertainties inherent to this approach. We present a detailed analysis of the impact of such uncertainties for four self-consistent axisymmetric CCSN models initiated from stellar metallicity, non-rotating progenitors of 12 $M_\odot$, 15 $M_\odot$, 20 $M_\odot$, and 25 $M_\odot$ and evolved with the smaller $\alpha$-network to more than 1 s after the launch of an explosion.

Tables

Table 1. Table 1: Initial distributions of Lagrangian tracer particles

	Models
Initial distribution	B12-WH07	B15-WH07	B20-WH07	B25-WH07
Number of particles ($N_{\mathrm{tp}}$)	4000	5000	6000	8000
Number of mass shells (rows)	100	125	150	200
Particles per mass shell (columns)	40	40	40	40
Mass per shell [$M_{☉}$$\times {10}^{-3}$]	7.472	11.46	14.18	13.94
Mass per particle ($\Delta m$) [$M_{☉}$$\times {10}^{-4}$]	1.868	2.864	3.545	3.486
Inner boundary
Inner row radius [km]	890.8	1,131	1,385	1,472
Inner edge mass [$M_{☉}$]	1.203	1.321	1.424	1.480
Average mean mass number ($\overline{A}$)	37.77	42.42	43.05	42.17
Outer boundary
Outer row radius [km]	15,000	19,456	19,751	19,751
Outer edge mass [$M_{☉}$]	3.879	5.476	3.551	4.269
Average mean mass number ($\overline{A}$)	15.16	15.01	16.99	17.05

Table 2. Table 2: Nucleosynthesis uncertainties

		Models
		B12-WH07	B15-WH07	B20-WH07	B25-WH07
$t_{\mathrm{pb}}$ at simulation end ($t_f$) [s]		1.336	1.200	1.383	1.399
Masses at $t_f$
Off-grid [$M_{☉}$]		8.829	10.04	12.32	11.26
Unshocked (on grid) [$M_{☉}$]		0.1900	0.3600	0.7852	1.368
$M_{\mathrm{unb}}$ [$M_{☉}$]		0.3964	0.7117	0.9919	1.284
$M_{+}$ [$M_{☉}$]		0.3893	0.6937	0.9334	1.052
$M_{-}$ [$M_{☉}$$\times {10}^{-2}$]		0.7659	1.804	5.849	23.14
$\left[M_{-}\right]$ [$M_{☉}$$\times {10}^{-2}$]		1.139	3.208	6.594	24.71
Combined uncertainties at $t_f$
$M_{-}/M_{\mathrm{unb}}$ (Section 4.1)		0.0193	0.0254	0.0590	0.1803
$\left[M_{-}\right]/M_{\mathrm{unb}}$ (Section 4.1)		0.0287	0.0451	0.0665	0.1926
$\Vert r_{{\delta }_{{\tau }^{\ast }}}\Vert$ (Section 4.2)		0.0043	0.0069	0.0101	0.0119
$\Vert r_{{\delta }_{\Delta 𝐫}}\Vert$ (Section 4.3)		0.6192	0.7117	0.7767	0.7086
$\Vert r_{{\delta }_{\mathrm{NSE}}}\Vert$ (Section 4.4)		0.1083	0.1134	0.0631	0.0770
Individual isotope uncertainties at $t_f$
$\left[{ϵ}^i\right]$ (sn150) (Section 4.1)	${}_{}{}^4\mathrm{He}$	0.1607	0.2070	0.3108	0.2678
	${}_{}{}^{28}\mathrm{Si}$	0.0863	0.1151	0.0502	0.3670
	${}_{}{}^{44}\mathrm{Ti}$	0.1722	0.1274	0.1956	0.1474
	${}_{}{}^{56}\mathrm{Ni}$	0.0701	0.0490	0.3239	0.3040
${\delta }_{{\tau }^{\ast }}^i$ (sn150) (Section 4.2)	${}_{}{}^{44}\mathrm{Ti}$	0.0272	0.0618	0.1000	0.0213
${\delta }_{\Delta 𝐫}^i$ ($\alpha$) (Section 4.3)	${}_{}{}^4\mathrm{He}$	0.0976	0.1762	0.1176	0.1600
	${}_{}{}^{28}\mathrm{Si}$	0.0248	0.0059	0.0117	0.0083
	${}_{}{}^{44}\mathrm{Ti}$	0.9383	0.9345	1.0703	1.0520
	${}_{}{}^{56}\mathrm{Ni}$	0.0121	0.0130	0.0247	0.0249
${\delta }_{\mathrm{NSE}}^i$ (sn150) (Section 4.4)	${}_{}{}^4\mathrm{He}$	0.2495	0.2682	0.3108	0.2236
	${}_{}{}^{28}\mathrm{Si}$	0.0000	0.0013	0.0502	0.0023
	${}_{}{}^{44}\mathrm{Ti}$	0.5531	0.5701	0.4101	0.4526
	${}_{}{}^{56}\mathrm{Ni}$	0.0350	0.0462	0.0437	0.0711

Equations

$$T_{\mathrm{NSE}}(\rho )=\{\begin{array}{ll}(\rho -{\rho }_{\mathrm{L}})\left(\frac{T_{\mathrm{H}}-T_{\mathrm{L}}}{{\rho }_{\mathrm{H}}-{\rho }_{\mathrm{L}}}\right)+T_{\mathrm{L}}& \text{if }\rho <{\rho }_{\mathrm{H}};\\ T_{\mathrm{H}}& \text{otherwise},\end{array}$$

$$r^{(n+1)}$$

$${\theta }^{(n+1)}$$

$$T(t)$$

$$\rho (t)$$

$$T(t)$$

$$\rho (t)$$

$$T(t)$$

$$\rho (t)$$

$$n({ϵ}_{\nu })=\frac{1}{F_2(\eta )T_{\nu }^3}\frac{{{ϵ}_{\nu }}^2}{\exp (({ϵ}_{\nu }/T_{\nu })-\eta )+1},$$

$$e_{\mathrm{tot}}\equiv e_{\mathrm{kin}}+e_{\mathrm{th}}+e_{\mathrm{grav}},$$

$$||r_{{\delta }_{{\tau }^{\ast }}}||\equiv \frac{{\sum }_i|{\delta }_{{\tau }^{\ast }}^i|}{{\sum }_i|{\log }_{10}(\sqrt{X_i({\tau }_{\max }^{\ast })X_i({\tau }_{\min }^{\ast })})|},$$

$$||r_{{\delta }_{\Delta \mathbf{r}}}||\equiv \frac{{\sum }_i|{\delta }_{\Delta \mathbf{r}}^i|}{{\sum }_i|{\log }_{10}(\sqrt{M_{\mathrm{C}}^i(t_{\mathrm{f}})M_{\mathrm{PP}}^i(t_{\mathrm{f}})})|}.$$

$$||r_{{\delta }_{\mathrm{NSE}}}||\equiv \frac{{\sum }_i|{\delta }_{\mathrm{NSE}}^i|}{{\sum }_i|{\log }_{10}(\sqrt{M_{\mathrm{PP}}^i(T_{\mathrm{NSE}}(\rho ))M_{\mathrm{PP}}^i(T_{\mathrm{NSE}}=8\mathrm{GK})})|}$$

Full text

methods: numerical — nuclear reactions, nucleosynthesis, abundances — stars: abundances — supernovae: general

††journal: ApJ

1 Introduction

The deaths of massive stars ($M>8\textrm{--}10$ $M_{\sun}$) as core-collapse supernovae (CCSNe) are an important link in our chain of origins from the Big Bang to the present. They are the dominant source of elements in the periodic table between oxygen and iron (Woosley & Weaver 1995; Thielemann et al. 1996), and there is strong evidence that they are correlated with the production of half the elements heavier than iron (Argast et al. 2004). Core-collapse supernovae serve both to disperse elements synthesized in massive stars during their lifetimes and to synthesize and disperse new elements themselves.

Modern multidimensional simulations by several groups, utilizing spectral neutrino transport, have successfully produced explosions for a variety of progenitors in axisymmetry (2D) though the explosion is delayed by hundreds of milliseconds compared to their non-spectral counterparts (Buras et al. 2003, 2006b, 2006a; Burrows et al. 2006, 2007; Bruenn et al. 2006, 2009; Marek & Janka 2009; Suwa et al. 2010; Takiwaki et al. 2012; Müller et al. 2012b, a; Bruenn et al. 2013; Müller & Janka 2014; Bruenn et al. 2016). Successful, fully self-consistent, spectral models have also achieved neutrino-driven explosions in three spatial dimensions (3D; Takiwaki et al. 2012, 2014; Melson et al. 2015a, b; Lentz et al. 2015). Qualitatively, the 2D and 3D simulations exhibit similar explosions, dominated by a small number of rising plumes that push the stalled shock outward. In both cases, explosions are inextricably linked to fluid instabilities that can only be properly accounted for in multidimensional simulations. Explosions can be aided by enhanced neutrino luminosities due to fluid instabilities in the proto-neutron star (proto-NS) (Smarr et al. 1981; Wilson & Mayle 1993), improved neutrino-heating efficiency behind the shock as a result of fluid motions induced by convective instabilities, and the standing accretion shock instability (SASI; Blondin et al. 2003; Blondin & Mezzacappa 2006; Foglizzo et al. 2007). However, 3D models exhibit a tendency to be delayed when compared to their 2D counterparts (Hanke et al. 2013; Melson et al. 2015b; Lentz et al. 2015; Janka et al. 2016). The exception is Melson et al. (2015b), wherein the authors evolve a qualitatively different progenitor (zero metallicity, 9.6 $M_{\sun}$) than the models discussed here and the works cited above.

Asymmetries are in fact present even in the pre-supernova core, with simulations going back two decades revealing complex, potentially interacting, convective carbon, neon, oxygen and silicon burning shells above the iron core (Bazan & Arnett 1998; Meakin & Arnett 2007; Arnett & Meakin 2011). Until recently, the effect of such asymmetries during core collapse, bounce and into the supernova explosion have been largely ignored. Couch & Ott (2013) and Müller & Janka (2015) demonstrated that artificially imposed but physically reasonable asphericities in the pre-collapse progenitor can qualitatively alter the evolution of the supernova after bounce. To build more realistic asphericities, Couch et al. (2015) subsequently completed the last few minutes of silicon shell burning in 3D before modeling the collapse and explosion of the star, also in 3D, with approximate neutrino transport. Further efforts in this direction seem warranted.

Despite the fundamentally multidimensional nature of a supernova explosion from its earliest moments, relatively limited work has addressed the impact of multidimensional behavior on the nucleosynthesis. The complexity and cost of self-consistent CCSN models (abetted by their past, frequent failures to produce explosions) has led the community to largely continue to rely on nucleosynthesis predictions for CCSNe from models using a parameterized kinetic energy piston (Rauscher et al. 2002; Woosley & Heger 2007; Magkotsios et al. 2010) or thermal energy bomb (Nagataki et al. 1998; Maeda et al. 2002; Umeda & Nomoto 2008). We have, however, learned that multidimensional effects produce significant differences in the fraction of ejecta which experiences $\alpha$-rich freeze-out (Nagataki et al. 1998; Maeda et al. 2002) and larger ejecta velocities, characterized by metal-rich clumps (Kitaura et al. 2006; Hammer et al. 2010; Ellinger et al. 2012; Wongwathanarat et al. 2013, 2015).

The work of Pruet et al. (2005), with its combination of a neutrino-driven explosion and multidimensional fluid flow, suggested that the impact on the nucleosynthesis of departing from stratified 1D simulations was significant. Nevertheless, very little investigation of CCSN nucleosynthesis from such models has been conducted, despite the multitude of exploding first-principles models with spectral neutrino transport in recent years. At least in part, this curious deficit can be attributed to the prolonged times after bounce the models must be evolved in order to fully characterize the ejecta and, therefore, compute the nucleosynthesis.

Nucleosynthesis studies of electron-capture supernovae (ECSNe), which arise from the collapse of oxygen-neon cores, are more mature, as these explosions trigger and complete more rapidly than in Fe-core SNe and can be obtained even in 1D simulations. Multi-dimensional investigations of ECSN nucleosynthesis (Wanajo et al. 2011, 2013a, 2013b) using spectral neutrino transport find only modest impact from multidimensional effects, which is unsurprising given the successful 1D explosions. Multi-dimensional Fe-core CCSN models also exhibit convective overturn near the outer proto-NS layers, potentially with even greater affect on the nucleosynthesis and, given the necessity of multidimensionality to engender these explosions, occurs merely as a subset of other multidimensional effects. Consequently, the lessons learned from ECSN nucleosynthesis studies, with respect to multidimensional effects, provide only modest insight into the CCSN problem.

This work aims to improve this insight by critically analyzing the sources of uncertainty in methods used to study CCSN nucleosynthesis in multidimensional models. In a subsequent paper, we will examine the detailed differences between the nucleosynthesis in these multidimensional models and parameterized, spherically symmetric (1D) models, but first we must verify that the uncertainties achieved in these multidimensional simulations are small enough that the comparison with 1D models is meaningful. To this end, in Section 2, we identify these uncertainties. In Section 3, we provide a description of the relevant numerical methods and key physical approximations employed by the our radiation-hydrodynamics code (Chimera) and subsequent post-processing calculations relating to the nucleosynthesis. In Section 4, we report our findings on the uncertainties at the end of our simulations, and in Section 5, we quantify how these have changed with time. Finally, in Section 6, we generalize our conclusions beyond the specifics of our models.

2 Nucleosynthesis Uncertainties

Calculating the nucleosynthesis of neutrino-driven, multidimensional CCSN models removes a range of necessary assumptions related to the various approaches used to generate explosions in parameterized investigations of this nucleosynthesis. These assumptions, though physically motivated, engender certain implicit uncertainties; foremost among these are the effects of neutrino-matter interactions and turbulent fluid flow. For fully self-consistent simulations, examination of the nucleosynthesis affords additional observable consequences of the explosion model that can be compared to observations—these may include, for example, detailed comparisons of the spatial distribution of radioactive nuclei (see, e.g., Wongwathanarat et al. 2016).

However, this progress to greater physical fidelity in our investigations of CCSN nucleosynthesis is not without compromise. Many of the bomb and piston models that have provided the bulk of our understanding of CCSN nucleosynthesis in recent years utilize realistic nuclear reaction networks and continue well after the supernova ejecta reaches the surface of the star, perhaps an hour after the formation of the proto-neutron star (see, e.g., Woosley & Heger 2007; Umeda & Nomoto 2008; Chieffi & Limongi 2013). To contain computational cost, multidimensional studies generally fail to include large nuclear reaction networks within their simulations. As a result, post-processing of Lagrangian thermodynamic histories with a realistic nuclear reaction network is required to generate abundances of all of the isotopes of interest. For Lagrangian methods, like smoothed particle hydrodynamics (SPH), these Lagrangian thermodynamic histories are simply the trajectories of the (active) particles that represent the fluid in this approach. For Eulerian methods, passive tracer particles must be added to the conventional hydrodynamic evolution, as these grid-based methods are otherwise unable to record the histories of fluid elements, instead evolving (and recording the histories of) fluid quantities at specified locations. Resolution is a limiting factor in any computational simulation. While SPH simulations, by their nature, provide Lagrangian fluid elements that fully sample the hydrodynamic evolution at the intrinsic resolution of the simulation, passive tracers within Eulerian simulations do not necessarily provide similar sampling.

Regardless of the formalism for the hydrodynamic evolution, the use of a realistic network only in post-processing removes the natural feedback between the evolution of the nuclear composition and the thermodynamic conditions. Thus, the accuracy of post-processing nucleosynthesis is only as good as the physical fidelity of the thermonuclear evolution included within the simulations. In most multidimensional CCSN studies, for matter at high temperature and density, this role is served by the nuclear equation of state under the assumption of nuclear statistical equilibrium (NSE), together with neutrino transport altering the neutronization. At lower temperatures, NSE does not apply; thus, a nuclear reaction network of some type is needed. The simplest approach used in such models is a flashing scheme, wherein “Oxygen” is converted directly to “Silicon” at some chosen temperature, and “Silicon” is similarly converted to NSE at another, higher temperature (see, e.g., Marek & Janka 2009). Of greater fidelity is an $\alpha$-network, composed of the 13 (or 14) isotopes with equal and even proton and neutron numbers between $\mathrm{{}^{4}He}$ and $\mathrm{{}^{56}Ni}$ (or $\mathrm{{}^{60}Zn}$) and the reactions, predominantly ($\alpha$,$\gamma$), that link them. While an $\alpha$-network is nearly complete for carbon burning and a reasonable approximation to oxygen burning, it does not include many of the reaction channels important for the production of iron and nickel (Timmes et al. 2000), either by silicon burning or by the recombination of dissociated nucleons and $\alpha$-particles, resulting in misestimation of the rate of nuclear energy release for these processes. A further limitation of an $\alpha$-network is the inability to track neutronization from electron and neutrino captures, so the neutronization must be monitored separately or ignored entirely. Even when handled separately, as it is in Chimera, it is difficult to evolve the neutronization correctly in the absence of NSE without all relevant species and reaction channels available.

The extreme cost of fully self-consistent simulations commonly limits our most realistic CCSN models to run only several hundred milliseconds after bounce, often aborting the simulations before the nucleosynthesis is complete. These cost-motivated compromises introduce additional sources of uncertainty in the calculation of CCSN nucleosynthesis that must be accounted for if we are to maximize the understanding gained from such studies. Early termination of these simulations requires extrapolation of the thermodynamic history if the thermonuclear evolution is still ongoing when the simulation stops. While a number of past studies have relied on similar extrapolations for their entire thermodynamic evolution, it is nonetheless important to constrain this uncertainty. Early termination also makes it more challenging to determine the mass-cut, which delineates the fates of the matter within the star, and judge the composition of the ejecta. The mass-cut does not impact all isotopes equally, the common perception (see, e.g. Diehl & Timmes 1998), based on spherically symmetric models, being that this uncertainty will most strongly impact the products of $\alpha$-rich freeze-out; however, this question merits re-examination in a multidimensional context.

3 Methodology

Defined by a shared set of numerical methods and physical approximations, the Chimera B-Series simulations originate from solar-metallicity, non-rotating progenitors calculated by Woosley & Heger (2007) for 12 $M_{\sun}$ (B12-WH07), 15 $M_{\sun}$ (B15-WH07), 20 $M_{\sun}$ (B20-WH07) and 25 $M_{\sun}$ (B25-WH07) stars. These simulations were evolved with the Chimera222http://www.ChimeraSN.org code described briefly herein, incorporating modern neutrino-matter interactions, self-consistent luminosities and neutrino spectra, and coupled nuclear burning. (For a more thorough discussion of methods in Chimera, see Bruenn et al. 2016). The B-Series simulations extend to 1.2–1.4 seconds after bounce and were carried out in axisymmetry from the onset of collapse, with the very small post-bounce roundoff errors supplying the perturbations that seed the growth of fluid instabilities. The long times after bounce to which the B-Series simulations were evolved make them a unique resource for investigations of the uncertainties stemming from ongoing hydrodynamic activity that can impact post-processing nucleosynthesis results, even at these late epochs. All times reported hereafter are relative to core-bounce, or equivalently as time post-bounce ($t_{\mathrm{pb}}$).

Chimera is a multidimensional, radiation-hydrodynamics code for stellar core collapse with five principal components: hydrodynamics, neutrino transport, self-gravity, a nuclear reaction network, and a nuclear equation of state (EoS). The details most relevant to this nucleosynthesis study are summarized below.

3.1 Radiation hydrodynamics

Hydrodynamics is evolved via a dimensionally split, Lagrangian-plus-remap scheme with piecewise parabolic reconstruction (PPMLR; Colella & Woodward 1984) as implemented in VH1 (Hawley et al. 2012), modified to include the consistent multi-fluid advection of Plewa & Müller (1999). The self-gravity is computed using a multi-pole expansion (Müller & Steinmetz 1995), replacing the Newtonian monopole with a general relativistic (GR) monopole (Marek et al. 2006, Case A). The neutrino transport solver is an improved and updated version of the multi-group (frequency), flux-limited diffusion (MGFLD) implementation of Bruenn (1985), which uses near-complete physics, solving for four neutrino species ($\nu_{e}$, $\bar{\nu}_{e}$, $\nu_{\mu\tau}=\{\nu_{\mu},\nu_{\tau}\}$, $\bar{\nu}_{\mu\tau}=\{\bar{\nu}_{\mu},\bar{\nu}_{\tau}\}$) while allowing for neutrino-neutrino scattering, pair exchange and other opacities.

3.2 Nucleosynthesis in Chimera

B-Series simulations utilize the $K=220$ MeV incompressibility version of the Lattimer & Swesty (1991) EoS for nuclear-matter densities ($\rho>10^{11}$ ${\mbox{g~{}cm}}^{-3}$) and an enhanced version of the Cooperstein (1985) EoS for $\rho<10^{11}$ ${\mbox{g~{}cm}}^{-3}$ where nuclear statistical equilibrium (NSE) applies. These provide a four-species representation of the chemical composition, consisting of free neutrons, free protons, $\alpha$-particles, and a representative heavy nucleus. To improve the fidelity of the composition of matter that may eventually become part of the ejecta, in regions of NSE where the electron fraction, $Y_{\mathrm{e}}$, is $\geq{}^{26}\!/_{56}$ (the value of $Z/A$ for $\mathrm{{}^{56}Fe}$), a 17-species representation of the composition is used, including free neutrons, free protons, the 14 even-Z and even-A nuclei between $\mathrm{{}^{4}He}$ and $\mathrm{{}^{60}Zn}$ which constitute an $\alpha$-network, as well as $\mathrm{{}^{56}Fe}$ to conserve $Y_{\mathrm{e}}$. Both the representative nucleus composition and Chimera’s 17-species NSE composition provide a fairly accurate representation of the nuclear composition, as long as NSE is maintained.

For regions not in NSE, the nucleosynthesis is computed within the constraints of an $\alpha$-network ($\alpha$, $\mathrm{{}^{12}C}$–$\mathrm{{}^{60}Zn}$) by XNet333http://eagle.phys.utk.edu/xnet/trac/, a fully implicit thermonuclear reaction network code (for details, see Hix & Thielemann 1999a). This network includes only one “effective” rate, triple $\alpha$, otherwise, only rates for reactions linking explicitly defined nuclei are used. For example, only ($\alpha$,$\gamma$) reaction rates are included in the $\alpha$-network, without the effective inclusion of ($\alpha$,p) and ($\alpha$,n) reaction rates. In regions evolved by the reaction network from the beginning of the simulations, the initial composition is determined by mapping the abundances of free nucleons and the $\alpha$-network species directly from the stellar progenitor to the computational grid in Chimera. An $\alpha$-network is of course quite limited in its ability to follow the evolution of composition where there should be a significant population of species with $N\neq Z$. To account for neutron-rich nuclei in the stellar progenitor, abundances of all heavy nuclei which are not part of the $\alpha$-network are bundled into an inert, representative auxiliary nucleus to conserve $Y_{\mathrm{e}}$. To account for free nucleons, the abundances of free nucleons are maintained in the non-NSE region, however they do not evolve.

Naturally, zones will transition in to and out of NSE during the simulation, as matter infalls toward the core, encounters the supernova shock and is expelled. The decision on when to transition from the NSE composition evolution to the nuclear reaction network, or vice versa, depends on trade-offs in the the ability to accurately calculate the composition and the computational cost. To determine whether a zone should transfer into NSE, and therefore be omitted from nuclear burning, Chimera applies an empirically determined linear relationship between the NSE transition temperature, $T_{\mathrm{NSE}}$, and density:

[TABLE]

where $\rho_{\mathrm{L}}\equiv 5\times 10^{7}~{}{\mbox{g~{}cm}}^{-3}$, $\rho_{\mathrm{H}}\equiv 2\times 10^{8}~{}{\mbox{g~{}cm}}^{-3}$, $T_{\mathrm{L}}\equiv 5.7~{}\mathrm{GK}$, and $T_{\mathrm{H}}\equiv 6.5~{}\mathrm{GK}$. This relationship is based on the need for NSE to be established, or maintained, over less than a Chimera timestep, otherwise, the assumption of NSE is in error. The expression for $T_{\mathrm{NSE}}(\rho)$ is based loosely on calculations of the peak temperatures required to achieve NSE in parameterized explosion simulations, starting with the one zone models of Woosley et al. (1973) and including 1D models (see, e.g., Thielemann et al. 1986, 1990). Operationally, at the beginning of a global timestep, any non-NSE zone for which $T\geq T_{\mathrm{NSE}}$ is transitioned to NSE, and the composition determined by the current temperature, density and $Y_{\mathrm{e}}$. Similar approaches are employed by many core collapse supernova simulations codes, though the expression for $T_{\mathrm{NSE}}$and the means to evolve the nuclear composition before NSE is achieved differ considerably.

The transition out of NSE is more challenging, as an initial composition must be generated for the nuclear reaction network. In fact, we find that the details of this transition out of NSE have a significant impact on the final composition of the ejecta (see Section 4.4). In Chimera, there are two cases, based on the ability to generate and evolve an accurate composition in the $\alpha$-network. These are formally delineated by whether the abundance of heavy nuclei would be dominated by $\mathrm{{}^{56}Ni}$, which is actively evolved by the network, or $\mathrm{{}^{56}Fe}$, which is held constant. If the composition would include more $\mathrm{{}^{56}Ni}$ than $\mathrm{{}^{56}Fe}$, effectively a $Y_{\mathrm{e}}\gtrsim 0.48$, a zone which is in NSE at the beginning of a timestep will be transitioned out of NSE if $T<T_{\mathrm{NSE}}-0.2~{}\mathrm{GK}$. In this case, the initial composition is generated by Chimera’s 17 species NSE calculation, fully populating the composition evolved by the $\alpha$-network, plus inert abundances of neutrons, protons and $\mathrm{{}^{56}Fe}$. From this detailed composition, continued compositional evolution, for example, the recombination of $\alpha$-particles into heavy nuclei, is possible. This is a key way in which Chimera’s treatment of the transition out of NSE is an improvement over most other core collapse supernova simulations codes, which have only the 4 species composition provided by the EoS with which to seed the network. For more neutron-rich conditions, the transition out of NSE does not occur until $T<4.9~{}\mathrm{GK}$. While this later transition runs the risk of suppressing the $\alpha$-rich freezeout under these conditions, the poor reproduction of neutron-rich compositions by the $\alpha$-network results in large, frozen populations of nucleons and neutron-rich heavy elements that are also unphysical. For this reason, and because lower $Y_{\mathrm{e}}$ is generally associated with higher densities and hence normal (or at least only mildly $\alpha$-rich) freezeout, the NSE composition provided by either Chimera’s 17 species NSE (for $Y_{\mathrm{e}}\geq{}^{26}\!/_{56}$) or the nuclear EOS is generally less wrong than that afforded by Chimera’s included network. Of course, the real solution to this problem is the use of a nuclear reaction network that is significantly more complete than the $\alpha$-network, a task we are undertaking (Lentz et al. 2017).

3.3 Tracer particle method

We provide a brief description of the implementation of The tracer particle method (TPM) in Chimera is implemented using passive, Lagrangian tracer particles to record the thermodynamic histories of individual mass elements that are then be used for post-processing nucleosynthesis calculations (see, e.g., Nagataki et al. 1997; Nishimura et al. 2006; Seitenzahl et al. 2010; Nishimura et al. 2015). Following each hydrodynamic directional sweep, the position of a tracer particle in that direction at time $t^{(n)}$, $(r^{(n)},\theta^{(n)})$, is advanced to $t^{(n+1)}$ according the simple Euler method, assuming constant velocity $(v_{r}^{(n)},v_{\theta}^{(n)})$ through the time interval $\Delta t^{(n)}=t^{(n+1)}-t^{(n)}$:

[TABLE]

Physical quantities are linearly interpolated in radius to tracer particle positions from the zone-center (cell-averaged) values of the computational grid; the lone exception being the interpolation of differential neutrino number fluxes, which are defined at radial zone edges.

The tracer particles are initially distributed into rows, radial shells uniformly-spaced in mass, beginning 0.1 $M_{\sun}$ inside the edge of the progenitor’s iron core. Thus, no tracers are placed in the inner iron core, conserving the tracers for matter more likely to be ejected. The particles within each row are placed into columns in latitude that represent uniform volume (and hence uniform mass): $\Delta(\cos\theta)=2/N$, where $N=40$ particles per row for B-Series models. For B12-WH07, $N_{\mathrm{tp}}=4,000$ tracer particles are used, translating to 40 columns and 100 rows. Each tracer represents $\Delta m=1.87\times 10^{-4}~{}M_{\sun}$, initially extending from $\approx$890 km to the carbon-enriched oxygen-shell at $\approx$15,000 km. This degree of tracer resolution is common among prior nucleosynthesis studies employing the TPM (Nagataki et al. 1997; Pruet et al. 2005; Nishimura et al. 2006) and is similar to the “medium” resolution case of Nishimura et al. (2015). For more massive progenitors, more tracers are employed. However, the number of tracers grows more slowly than the mass in the silicon and oxygen shells, resulting in lower mass resolution in the more massive models. A description of the initial tracer particle distribution and representative mass of each tracer particle for each of the B-Series models is given in detail in Table 1.

3.3.1 Parameterized thermodynamic trajectories

Purely parameterized expansion profiles have long been used to study CCSN nucleosynthesis. Fowler & Hoyle (1964) introduced the use of adiabatic expansion following a peak temperature, $T_{\mathrm{peak}}$, and peak density, $\rho_{\mathrm{peak}}$, mimicking the passage of the shock-wave, with a characteristic expansion timescale, $\tau^{*}$, equal to the free-fall timescale ($446/\sqrt{\rho_{\mathrm{peak}}}$ s if $\rho_{\mathrm{peak}}$ is measured in ${\mbox{g~{}cm}}^{-3}$). With this assumption, the thermodynamic trajectory for a single particle may be expressed as

[TABLE]

where $\Delta t=t-t_{\mathrm{peak}}$ is the time since the shock passage.

A number of explorations (see, e.g, Woosley et al. 1973; Meyer et al. 1998; Hix & Thielemann 1999b) have used this simple parametrized model for their entire thermodynamic evolution. Magkotsios et al. (2010) utilized both this exponential decline as well as the power-law decline that results from the homologous expansion of a uniform sphere. Nishimura et al. (2015) also employ a power-law description of the expanding matter for simulations of magnetically-driven, rotating CCSNe soon after the early ejecta drops below NSE conditions. Panov & Janka (2009) employed an exponential decrease in temperature and density to describe the first (adiabatic) stage of a homologously expanding neutrino-driven wind before switching to a power-law model at later times, representative of reduced wind acceleration occurring after early homologous expansion. In this case, the density and temperature decline much less steeply during the power-law phase than during exponential expansion and can be calculated from the density ($\rho_{0}$) and temperature ($T_{0}$) at the time $t_{0}$ when the matter is subjected to the reduced acceleration. For $t>t_{0}$,

[TABLE]

Despite their efforts to account for different phases in the expansion, these purely analytic parameterizations usually begin at peak temperatures and densities and expand smoothly, failing to directly capture thermodynamic variations seen in multidimensional models.

Parameterizations like these can also be employed to extrapolate self-consistent thermodynamic trajectories from hydrodynamic simulations to much later times than the original simulations. In that case, peak values of temperature and density are replaced, for example in Equation 5, by final values, $T_{\mathrm{f}}$ and $\rho_{\mathrm{f}}$, and $\Delta t$ becomes $t-t_{\mathrm{f}}$, the time beyond the end of the simulation. With the preceding hydrodynamic evolution in hand, better estimated expansion timescales are possible. For the particle data used herein, the effective expansion timescale, $\tau^{*}$, is calculated by averaging the instantaneous value, $\tau=-\rho/\dot{\rho}$, during periods of expansion near the end of the simulation. The thermodynamic histories of particles are often quite noisy, so a 25-ms wide, six-degree Savitzky-Golay smoothing filter is applied to the density profile for computation of numerical derivatives (Press et al. 2007). For the epochs we discuss in this paper, we find that the final 150 ms of the simulation provides an adequate time window for properly sampling the particle behavior. The time window is also narrowed such that deviations in $T^{3}/\rho$ from the final value do not exceed 10%. Furthermore, we exclude outliers in the calculation of $\tau^{*}$ by only considering the interquartile range of discrete values of $\tau$ over the specified time window.

The exponential behavior in temperature and density is supplemented with a more slowly declining power-law trajectory (Equation 7) after nuclear burning has reached freeze-out conditions ($T\lesssim 0.5$ GK) at time $t_{\mathrm{fo}}$. This prevents the small temperatures that result from the exponential expansion from producing spurious results below the REACLIB lower temperature limit. The extrapolated thermodynamic trajectory for each tracer particle combines both of these parameterizations:

[TABLE]

For neutrino-induced reactions extrapolated to $t>t_{\mathrm{f}}$, we assume that $T_{\mathrm{\nu}}$ remains constant and the neutrino luminosity is constant therefore $\phi_{\mathrm{\nu}}\propto 1/r(t)^{2}$, where $r(t)=r(t_{\mathrm{f}})+v_{\mathrm{r}}(t_{\mathrm{f}})(t-t_{\mathrm{f}})$ assuming constant radial velocity.

3.4 Nuclear reaction network

Beyond the $\alpha$-network described in Section 3.2, we also utilize a 150-species network, which we will refer to as the sn150-network or just sn150. Ranging in atomic number from $Z=0$ to $Z=30$ and including all isotopes from neutron number $N=Z$ to the most neutron-rich stable isotope for each element, sn150 represents a first-order improvement to the $\alpha$-network. While not sufficient to capture more exotic nuclear processes like the ${\nu}p$-process and $r$-process, this moderately-sized reaction network encompasses a significant fraction of elemental abundances and energy-producing reactions important to the core-collapse problem, allowing proper neutronization and, when coupled to the hydrodynamics, a more accurate nuclear energy generation rate. Reaction rates are taken from the REACLIB444https://groups.nscl.msu.edu/jina/reaclib/db/ compilation (Cyburt et al. 2010) and supplemented/supplanted with $\beta$-decay rates and electron capture rates on free nucleons and heavy nuclei (Fuller et al. 1985; Oda et al. 1994; Langanke & Martínez-Pinedo 2000).

Neutrino capture rates on free nucleons and nuclei are also included in post-processing calculations. For this purpose, the integrated number flux of the neutrino distribution, $\phi_{\mathrm{\nu}}$, is recorded for $\nu_{e}$ and $\bar{\nu}_{e}$ at the tracer particles’ location. Fluxes for $\nu_{\mu\tau}$ and $\bar{\nu}_{\mu\tau}$ are also recorded but are not included in the post-processing nucleosynthesis. The neutrino-induced reaction cross-sections also require the neutrino temperatures, $T_{\mathrm{\nu}}$, which we calculate by fitting the numerical neutrino distribution from the radiation-hydrodynamics simulation, $n(\epsilon_{\nu})$, to a Fermi-Dirac spectrum of arbitrary degeneracy $\eta$:

[TABLE]

where $F_{2}(\eta)$, the second-order Fermi-Dirac integral, is used for normalization to unity.

Post-processing calculations for ejected tracer particles which attain NSE begin at a point in time near the transition out of NSE ($t_{\mathrm{NSE}}$), with initial conditions determined by finding the appropriate NSE composition which solves the Saha equation for the selected network and the corresponding thermodynamic conditions: $\rho(t_{\mathrm{NSE}})$, $T(t_{\mathrm{NSE}})$, and $Y_{\mathrm{e}}(t_{\mathrm{NSE}})$ (see, e.g., Hix & Meyer 2006; Hartmann et al. 1985, and references therein, for more details on the NSE equations and resulting abundances). For tracers that never reach NSE, the initial abundances, temperature, and density are interpolated directly from the composition provided in the stellar progenitor. We then evolve the composition along the TPM-generated thermodynamic trajectories with a standalone implementation of XNet using a fully implicit integration scheme (Hix & Thielemann 1999a). In this way, we are able to generate more detailed nucleosynthetic yields using more realistic nuclear reaction networks.

4 Quantifying Uncertainties

We quantify contributions to uncertainty in multidimensional models of CCSN nucleosynthesis relating to the determination of the multidimensional mass-cut (Section 4.1), estimation of expansion timescales for use in thermodynamic extrapolation (Section 4.2), and limitations in the spatial (Section 4.3) and temporal (Section 4.5) resolution of the Lagrangian tracer particles. In Section 4.4 we discuss the impact that choices in the transition from NSE has on the composition. We also examine some consequences of the nuclear reaction network size. The full impact of using an $\alpha$-network must await similar simulations with a large in situ nuclear reaction network. Some additional sources of uncertainty, for example, those related to the stochasticity inherent in multidimensional models, are beyond the scope of this study.

4.1 Determination of the multidimensional mass-cut

Observations and models agree that CCSNe are highly asymmetric events driven by complex and/or turbulent fluid flows. The implications of this inherent multidimensionality on the nucleosynthesis are lost in 1D simulations, wherein a clear distinction, the mass-cut, is easily made between matter which is ejected to the interstellar medium and that which falls back to the proto-NS (defined hereafter as the region where $\rho>10^{11}$ ${\mbox{g~{}cm}}^{-3}$). While the placement of the mass-cut can vary in time in a 1D model, for example, as fallback results from matter formerly expanding being decelerated by interactions with the stellar envelope, the mass-cut is nonetheless a unique mass coordinate within the progenitor star. Extending this distinction to 2D and 3D simulations is challenging, as there is no requirement that the ejecta form a contiguous region, an implicit feature of 1D simulations. Further, the ability for accretion and outflow to occur simultaneously makes the determination of even the initial mass-cut a much more gradual process.

Ultimately, this requires evolving a model well beyond the initial development of an explosion until such time that the downflows, which have long been cut-off from the rest of the star at the shock, cease falling onto the proto-NS and accretion dramatically slows. Here, the limitation in our ability to extend the models, due to the computational cost of running simulations with spectral neutrino transport, is a great impediment compared to typical parameterized models. In this study, we have extended four self-consistent models with spectral neutrino transport using Chimera much further than similar models have heretofore been run in order to examine this issue. We find that the initial multidimensional mass-cut is not yet fully resolved, despite 1.2–1.4 s of evolution after bounce and $\approx$1 s after the initiation of an explosion.

Examining the mass-cut begins by defining the ejecta. Once the supernova shock breaks through the surface of the star, the ejecta is self-defined as the matter which propagates into the ISM. However, a definition is needed that can be applied at much earlier epochs. The treatment of the explosion energy faces a similar challenge of definition at these early times; thus, we mirror here the treatment of the explosion energy discussed in Bruenn et al. (2016), which is similar to the approach used by several authors in recent years. With contributions from specific kinetic ($e_{\mathrm{kin}}$), thermal ($e_{\mathrm{th}}$), and gravitational binding energy ($e_{\mathrm{grav}}$), the specific total energy,

[TABLE]

defines the unbound ejecta as particles for which $e_{\mathrm{tot}}>0$. We label this set of unbound particles as $\mathbb{P}_{\mathrm{unb}}$ and the corresponding mass as $M_{\mathrm{unb}}$. This should not to be confused with the total ejected mass, $M_{\mathrm{ej}}$, which also includes contributions from as yet unshocked matter, both on and above Chimera’s computational grid. Likewise, we use $\mathbb{P}_{\mathrm{bound}}$ and $M_{\mathrm{bound}}$ to represent the bound matter ($e_{\mathrm{tot}}<0$) outside of the proto-NS. Ideally, all tracer particles in $\mathbb{P}_{\mathrm{unb}}$ would constitute a portion of the observed ejecta. However, due to the work required to lift the stellar envelope out of the star’s gravitational well, $\mathbb{P}_{\mathrm{unb}}$ is likely an overestimate of particles that would ultimately be ejected.

The free-fall time, $t_{\mathrm{ff}}=\sqrt{3\pi/(32G\rho)}$ s (Kippenhahn & Weigert 1990), provides a rough estimate of a lower-bound for the additional time needed to establish the initial multidimensional mass-cut (i.e. $\delta t_{\mathrm{masscut}}\gtrsim t_{\mathrm{ff}}$). In addition to the free-fall time, each parcel of matter represented by a particle is subject to additional neutrino heating as it nears the proto-NS and may be re-ejected. Letting $\rho$ be the average density of all marginally unbound ($|e_{\mathrm{grav}}|\leq e_{\mathrm{kin}}+e_{\mathrm{th}}\leq 1.5|e_{\mathrm{grav}}|$) particles with relatively small radial velocities ($v_{\mathrm{r}}<1\times 10^{3}~{}\mbox{km~{}s}^{-1}$), we estimate lower-bounds for $t_{\mathrm{f}}+\delta t_{\mathrm{masscut}}$ as 6.50 s, 4.29 s, 4.37 s, and 4.32 s for B12-WH07, B15-WH07, B20-WH07, and B25-WH07, respectively.

A further complication is that some particles in $\mathbb{P}_{\mathrm{unb}}$ have negative radial velocities and, thus, increasing temperature and pressure. Therefore, they cannot be reliably extrapolated for post-processing nucleosynthesis. For this reason, it is helpful to categorize $\mathbb{P}_{\mathrm{unb}}$ by radial velocity, which we label as $\mathbb{P}_{+}$ and $\mathbb{P}_{-}$ for $v_{\mathrm{r}}>0$ and $v_{\mathrm{r}}<0$, respectively. We find that the ultimate fates of particles in $\mathbb{P}_{-}$ as either ejecta or part of the proto-NS ($\mathbb{P}_{\mathrm{PNS}}$) are often unknowable at the end of the simulation. Consequently, the mass represented by $\mathbb{P}_{-}$, $M_{-}$, is one indication of uncertainty in the total ejecta mass. Therefore, trends in $M_{-}(t)$ partially characterize the duration of time we must evolve a model in order to keep this uncertainty manageable as illustrated in Figure 1 (red lines). Total unbound mass ($M_{\mathrm{unb}}$, black lines) gradually increases as the shock lifts the envelope. $M_{\mathrm{bound}}$ (green lines) remains stubbornly above 0.001 $M_{\sun}$ at the end of each model, reflecting the ongoing accretion evident in these models. $M_{-}$ (blue lines) still exhibits significant growth at the end of B20-WH07 and B25-WH07, while B12-WH07 and, possibly, B15-WH07 have leveled out near 0.01 $M_{\sun}$. We can attribute this behavior to an immature explosion in B25-WH07 and a weak explosion in B20-WH07. Across all four models, the final values of $M_{-}$ exhibit a monotonic increase with progenitor mass from 0.008 $M_{\sun}$ in B12-WH07 to 0.2 $M_{\sun}$ in B25-WH07.

For B12-WH07, $M_{-}$ is relatively unchanged in the 100 ms prior to the end of the simulation, defined as $t_{\mathrm{f-100}}$. However, closer inspection reveals that, while the total number of particles in $\mathbb{P}_{-}$ tracks the behavior of $M_{-}$, the individual particles in $\mathbb{P}_{-}$ are changing as convective flows and shear move particles between $\mathbb{P}_{-}$ and $\mathbb{P}_{+}$. This variability of $v_{\mathrm{r}}$ for $\mathbb{P}_{\mathrm{unb}}$ near the end of the simulation complicates the determination of a particle as being representative of the ejected matter. In Figure 2, we illustrate this phenomenon by plotting the various fates of the particles in $\mathbb{P}_{-}(t)$ (black line) at $t-100$ ms.

Ideally, $\mathbb{P}_{-}(t)$ and $\mathbb{P}_{-}(t-100~{}\mathrm{ms})$ would be identical, indicating a consistent determination of the particles’ fates. However, since there is a persistent fraction of particles in $\mathbb{P}_{-}(t)$ that are being classified as ejecta only 100 ms in the past ($\mathbb{P}_{+}(t-100~{}\mathrm{ms})$; red line), we identify all particles that were in $\mathbb{P}_{-}$ at any time in the last 100 ms of the simulation, which we label $[{\mathbb{P}_{-}}]$, as having indeterminate fates and contributing to a better estimated uncertainty in the ejecta mass. We use this type of nomenclature hereafter to refer to quantities derived from $[{\mathbb{P}_{-}}]$ (e.g. $[{M_{-}}]$). For B12-WH07, we see that $[{M_{-}}]$ is $\sim$50% larger than $M_{-}$, reflecting that a significant portion of $\mathbb{P}_{-}$ is undergoing substantial variations in radial velocity over the past 100 ms. For B15-WH07, $[{M_{-}}]/M_{-}=1.78$ is relatively even larger, indicating more exchange of inwardly and outwardly moving particles. Intriguingly, the much larger and still growing number of particles in $\mathbb{P}_{-}$ in B20-WH07 and B25-WH07 are not undergoing similar variations, with $[{M_{-}}]$ only 13% and 7% larger than $M_{-}$, respectively.

The fates and peak temperatures, $T_{\mathrm{peak}}$, for shocked particles in the inner 8,000 km of each of the B-Series models are shown at $t_{\mathrm{f}}$ in Figure 3. The shock has generally moved beyond this grid, except in the upper left corners for B20-WH07 and B25-WH07, wherein the absence of shocked tracers is visible. Downflows of bound matter ($\mathbb{P}_{\mathrm{bound}}$; open, gray circles) are clearly visible, especially in B12-WH07 and B15-WH07, dragging unbound but now infalling particles ($\mathbb{P}_{-}$; open, colored circles). The much larger mass represented by $\mathbb{P}_{-}$ particles, $M_{-}$, for B20-WH07 and B25-WH07, illustrated in Figure 1, is plainly visible in the open, colored circles of Figure 3, especially for B25-WH07. Clearly, in both B20-WH07 and B25-WH07, the shock in the equatorial regions has been less successful in pushing the envelope outward, suggesting significant accretion will continue for a considerable time.

Since the tracers are a Lagrangian representation of the matter, we can trace them back in time to their original positions. Figure 4 shows the same quantities as Figure 3, but with the tracers positioned at their initial locations, revealing the fate of tracers as a function of their starting point in the progenitor. These same progenitor models were exploded by Woosley & Heger (2007) using a parameterized piston positioned at the $S/N_{A}\mathrm{k_{B}}=4.0$ isoentropy contour. This represents the mass-cut in these parameterized models, $M_{\mathrm{cut}}^{\mathrm{1D}}$, and generally lies at the inner boundary of the oxygen-burning shell (approximately 2,800 km, 3,900 km, 2,600 km, and 2,800 km for the 12 $M_{\sun}$, 15 $M_{\sun}$, 20 $M_{\sun}$, and 25 $M_{\sun}$ B-Series models, respectively). The location of $M_{\mathrm{cut}}^{\mathrm{1D}}$ typically lies farther from the iron-core than its 2D counterpart, $M_{\mathrm{cut}}^{\mathrm{2D}}$, calculated as the combined mass of proto-NS material and bound material at $t_{\mathrm{f}}$. It is important to note, however, that while a spherical mass-cut can be constructed for 2D models, it is not a good representation of the fate of individual mass elements. While detailed analysis of the differences between spherical, parameterized models and multidimensional, neutrino-driven models will be presented in Harris et al. (2017), the presence of bound particles originating above the spherical mass-cut, and unbound particles below, illustrates the fundamental difference in a multidimensional mass-cut. Intriguingly, a significant number of ejected tracers originate in the relatively dense silicon shell, some nearly to the edge of the iron core, which has ramifications for the nucleosynthesis. Of more immediate interest is the identification of the currently accreting matter.

In B12-WH07, B15-WH07, and B20-WH07, $\mathbb{P}_{-}$ (open, colored circles) covers a small, isolated region of the initial particle distribution (see Figure 4), though there is a trend to wider regions in latitude with increasing mass. These particles represent an uncertainty in the identification of the ejecta (i.e. mass-cut) and originate 3,000–4,000 km from the center and have reached peak temperatures of 3–5 GK. From this, we can infer that the explosive burning products of a composition originally consisting of $\mathrm{{}^{16}O}$ and $\mathrm{{}^{28}Si}$ are most susceptible to the mass-cut. This is in contrast to the usual assumption that products of $\alpha$-rich freeze-out are most susceptible to this determination. The characteristics of $\mathbb{P}_{-}$ in B25-WH07 are much more extreme than that of the other three models. In the case of the 25 $M_{\sun}$ model, the uncertainty in the identification of the ejecta is not confined to a small region in the initial particle distribution, but nearly extends from pole-to-pole in latitude.

We quantify these uncertainties in all four models in Figure 5 by comparing the total unbound mass of each isotope after post-processing with sn150 to that represented by $\mathbb{P}_{-}$ and $[{\mathbb{P}_{-}}]$. The values of $\epsilon\equiv M_{-}/M_{\mathrm{unb}}$ and $[{\epsilon}]\equiv[{M_{-}}]/M_{\mathrm{unb}}$ for each model, as well as $[{\epsilon^{i}}]$ for $\mathrm{{}^{4}He}$, $\mathrm{{}^{28}Si}$, $\mathrm{{}^{44}Ti}$, and $\mathrm{{}^{56}Ni}$, are given in Table 2. The larger disagreement of $\epsilon$ relative to $[{\epsilon}]$ in B12-WH07 and B15-WH07 compared to the less evolved B20-WH07 and B25-WH07 models relates to the different behavior of $M_{-}$ and $[{M_{-}}]$ exhibited in B20-WH07 and B25-WH07, and is a testament to the challenge of running simulations sufficiently past bounce such that the multidimensional mass-cut is truly resolved.

As illustrated by Figure 5, the effect of the indeterminate mass-cut in B12-WH07 is most prevalent for $A\geq 28$, and $[{\epsilon^{i}}]\lesssim 0.2$. The impact on the production of these species can be understood by considering the region of the star where the mass-cut has not yet been determined, illustrated by the open, colored circles in Figure 3. For B12-WH07, this region is confined to the inner 5,000 km of the star around a cut-off downflow rich in $\mathrm{{}^{28}Si}$ which continues to accrete onto the proto-NS long after the development of the explosion. This region is similarly confined for both B15-WH07 and B20-WH07, though characterized by more matter which has had its initial composition effectively reset by attaining NSE, evidenced by greater uncertainties for iron-group isotopes. This trend also appears in Figure 5, where $\epsilon^{i}$ in B15-WH07, B20-WH07, and B25-WH07 is characteristically different from B12-WH07 for $A>48$, suggesting that the isotopes of Cr, Mn, Fe, and Co with values of $Z/A$ closest to $Y_{\mathrm{e}}(t_{\mathrm{NSE}})$ of the ejected material are especially susceptible to the partially-resolved, multidimensional mass-cut. In Figure 3, $\mathbb{P}_{-}$ in B25-WH07 spans a much larger region of the star, not only the result of a deeper gravitational well, but also a consequence of an explosion that has been less successful in lifting the marginally-unbound, equatorial downflow (see animated Figure 2 in Bruenn et al. 2016).

4.2 Thermodynamic extrapolation

The distribution of particle temperatures at $t_{\mathrm{f}}$ in each B-Series simulation, $T_{\mathrm{f}}$, is shown in Figure 6 (red) for $\mathbb{P}_{+}(t_{\mathrm{f}})$. In the case of the relatively less-evolved B25-WH07 (see Figure 12 in Bruenn et al. 2016), $\approx$0.11 $M_{\sun}$ of the ejecta continues to experience explosive burning ($T_{\mathrm{f}}\gtrsim 3$ GK), despite 1.4 s of evolution after core-bounce. For the other three models, the nuclear reactions which account for the bulk of the nickel production in CCSN ejecta have ceased. However, secondary nuclear burning processes will continue to alter the abundance distribution until the matter freezes out (Woosley et al. 1994), and proton and neutron captures will continue until the temperature of the ejecta falls below $\approx$0.5 GK (Fröhlich et al. 2006), particularly in proton-rich ejecta. Short of being able to extend the simulations to this freeze-out temperature, we must extrapolate the thermodynamic conditions using the method described in Section 3.3.1.

In order to determine the collective effect of single-particle extrapolation uncertainties, we post-process thermodynamic profiles generated by extrapolating from points in time in the last 150 ms of the simulation corresponding to the minimum and maximum estimates of the expansion timescales, $\tau^{*}_{\min}$ and $\tau^{*}_{\max}$, respectively. These changing timescales reflect the effects of ongoing hydrodynamic activity. Of course, any extrapolation will fail to capture future hydrodynamic activity that deviates from true isentropic expansion, but changes over the last 150 ms are a reasonable proxy to estimate this uncertainty. Defining a residual norm for the final composition,

[TABLE]

where $\delta_{\tau^{*}}^{i}=\log_{10}(X_{i}(\tau^{*}_{\max})/X_{i}(\tau^{*}_{\min}))$ and $X_{i}$ is the mass fraction of species $i$, we are able to easily identify particles whose nuclear products are particularly susceptible to such activity. Tracer particles with the largest extrapolation uncertainties, as measured by $||r_{\delta_{\tau^{*}}}||$, typically fall very close to the proto-NS surface and reach peak temperatures in excess of 20 GK before being ejected at high speed (see Figure 7). The exposure of such particles to large neutrino fluxes near the neutrino-emitting surface, the neutrino-sphere, is consistent with our premise that thermodynamic extrapolations initiated from $T_{\mathrm{f}}\lesssim 2$ GK are necessitated by the ongoing proton and neutron captures that typify neutrino-induced nucleosynthesis. In Figure 8, we show differing estimates of the expansion timescale for four particles from B12-WH07 with the largest values of $||r_{\delta_{\tau^{*}}}||$: B12-WH07-P1289, B12-WH07-P1422, B12-WH07-P1616, and B12-WH07-P1737. The differences are especially large for particle B12-WH07-P1289, which experiences a brief period of heating just prior to $t_{\mathrm{f}}$, pushing the temperature above 3 GK and reengaging nuclear reactions that had otherwise ceased. Hydrodynamically, this corresponds to the particle briefly halting its expansion as it retreats back towards the proto-NS before being swept up in the ejecta once more (see Figure 9).

However, this specific type of profile is very uncommon—particles with $||r_{\delta_{\tau^{*}}}||>0.1$ make up less than 1% of the ejecta mass in B12-WH07. In fact, B12-WH07-P1289 would have been a member of $\mathbb{P}_{-}$ just 50 ms earlier. For B12-WH07-P1422, B12-WH07-P1616, and B12-WH07-P1737, the differences in the thermodynamic extrapolations are less drastic but not insignificant. In each of these cases, the differing estimates of the expansion timescale can be attributed to hydrodynamical flows which nudge the particles in question, but not to the extent that expansion ceases outright (see Figure 9). In so doing, a brief period of heating breaks the assumption of isentropic expansion. For example, B12-WH07-P1422 exhibits significant differences in the final composition as a result of extrapolations that are quite similar but lagged by $\approx$150 ms. Estimates of the expansion timescale in B12-WH07-P1737 differ by a factor of four despite the extrapolations initiating only $\approx$50 ms apart. B12-WH07-P1616 represents a peculiar case; the alteration to its spatial trajectory, evident in Figure 9, sends it on a collision course with the persisting equatorial downflow. Though the particle is considered part of the ejecta at $t_{\mathrm{f}}$, it is now part of a convective eddy and will inevitably soon transition to $\mathbb{P}_{-}$. Not only does this highlight how nucleosynthesis predictions based on analytic extrapolations can fail to directly capture future hydrodynamic activity, but also elucidates the more subtle effect of different estimates of the expansion timescale.

Figure 9 illustrates the extent to which the spatial trajectories of these tracer particles can differ from those of their immediate neighbors, highlighting the highly asymmetric behavior which complicates multidimensional nucleosynthesis calculations. As previously mentioned, a shared yet rare characteristic of the selected particles (B12-WH07-P1289, B12-WH07-P1422, B12-WH07-P1616, and B12-WH07-P1737) is their exposure to large neutrino fluxes near the neutrino-sphere. As evidenced by their immediate neighbors, most tracer particles which descend deep into the gravitational well of the proto-NS become bound and do not represent ejecta matter.

At 100 s after bounce, the composition resulting from these varying extrapolations can be markedly different. As shown in Figure 10, deviations from predicted trajectories after $t_{\mathrm{f}}$ can lead to non-trivial uncertainties in the composition for individual particles.

The tracers discussed above were selected by their extreme uncertainties. For most tracers, and therefore the total ensemble, the impact of the uncertainties is much less pronounced. The extent to which this impacts the total ejecta mass for each nuclear species $i$, $M_{+}^{i}$, represented by $\mathbb{P}_{+}$ (the particles for which we perform extrapolations) is shown for each B-Series model in Figure 11. The relatively low temperatures of particles in $\mathbb{P}_{+}(t_{\mathrm{f}})$ (see Figure 6) helps to limit the overall impact of extrapolation uncertainties to only a few nuclear products. In general, $\delta_{\tau^{*}}^{i}\equiv\log_{10}(M_{+}^{i}(\tau^{*}_{\max})/M_{+}^{i}(\tau^{*}_{\min}))$ is larger for low-yield isotopes (e.g. $\mathrm{{}^{3}He}$, $\mathrm{{}^{13}C}$, $\mathrm{{}^{18}O}$, …), wherein the extrapolation uncertainties for a small subset of tracer particles can have a greater effect. We also identify trends of larger uncertainties for some isotopes that can’t be entirely attributed to a small number of tracer particles (e.g. $\mathrm{{}^{38}K}$, $\mathrm{{}^{43}Sc}$, $\mathrm{{}^{56}Co}$, and $\mathrm{{}^{57}Co}$). Uncertainties in $\mathrm{{}^{44}Ti}$ production, important to supernova remnant observations, are small, but non-negligible ($0.79\lesssim M_{+}^{i}(\tau^{*}_{\max})/M_{+}^{i}(\tau^{*}_{\min})\lesssim 0.95$; see Table 2 for values).

Other isotopes are affected by variations in $\tau^{*}$, but their identities vary between models. B20-WH07, for example, is uniquely characterized by a noticeable increase in $\delta_{\tau^{*}}^{i}$ for $12\leq A\leq 30$. B25-WH07 is also unique among the B-Series simulations in that the extrapolation error can be larger than 20% for several of the most neutron-rich isotopes of elements for which $Z\geq 14$, including $\mathrm{{}^{32}Si}$, $\mathrm{{}^{36}S}$, $\mathrm{{}^{40}Ar}$, $\mathrm{{}^{49}Sc}$, $\mathrm{{}^{50}Ti}$, $\mathrm{{}^{51}V}$, $\mathrm{{}^{53}Cr}$, $\mathrm{{}^{54}Mn}$, $\mathrm{{}^{58}Fe}$, and $\mathrm{{}^{59}Co}$.

As a single measure of extrapolation uncertainties for a given model, we define a global residual using Equation 12 but replacing the mass fraction, $X_{i}$, with $M_{+}^{i}$. Not surprisingly, $||r_{\delta_{\tau^{*}}}||$ increases with progenitor mass but is small relative to the other uncertainties (see Table 2).

4.3 Spatial resolution

Any simulation, by its nature, replaces continuous thermodynamic variables with a discrete number of elements, either in space (Eulerian) or in mass (Lagrangian). This discretization adds an uncertainty that is magnified if the number of elements is too small to capture the important features in sufficient detail. For post-processing nucleosynthesis, the resolution in question is the number and distribution of the Lagrangian tracer particles used to gather the thermodynamic histories.

The initial distribution of tracer particles in our simulations guarantees that we are sampling the ejecta uniformly in mass. Consequently, there are density-dependent variations in the spatial resolution, $\Delta\mathbf{r}$, of the ejecta. If sufficiently large, these variations may limit the ability of Lagrangian tracer particles to sample the ejecta at resolutions which reproduce the nuclear burning conditions initially encountered in the simulation. In particular, tracer particles tend to under-resolve regions of relatively low density. While we could adjust the initial tracer mass to better sample low-density regions in the progenitor, our ability to do this for dynamically developing low-density regions is limited. Past studies of the convergence properties of tracer particles in supernova nucleosynthesis (see, e.g., Seitenzahl et al. 2010) have been limited to Type Ia supernovae, which exhibit a narrower entropy range.

We introduce here a novel way to explore the limitations of tracer post-processing by post-processing with a network identical to that used in situ within Chimera. Care is taken to also generate initial abundances, both for particles that reach NSE and those that do not, that are identical to the methods used within Chimera. Thus, we have two representations of the nuclear composition that differ only in their effective resolution and the inability of the Lagrangian tracer particles to mix their composition. The total ejected mass of each nuclear species $i$ from in situ nucleosynthesis, $M_{\mathrm{C}}^{i}$, is calculated by integrating over all zones where $e_{\mathrm{tot}}>0$. The ejected mass from the equivalent post-processing calculation, $M_{\mathrm{PP}}^{i}$, is calculated with the $e_{\mathrm{tot}}>0$ criterion applied to particles instead of zones. The relative difference between in situ and post-processing calculations for individual isotopes, $\delta_{\Delta\mathbf{r}}^{i}\equiv\log_{10}(M_{\mathrm{PP}}^{i}/M_{\mathrm{C}}^{i})$, $M_{\mathrm{C}}^{i}(t)$, and $M_{\mathrm{PP}}^{i}(t)$ are shown throughout each B-Series simulation in Figure 12.

There is general agreement between $M_{\mathrm{C}}^{i}(t)$ and $M_{\mathrm{PP}}^{i}(t)$ for most species. $\mathrm{{}^{28}Si}$ and $\mathrm{{}^{56}Ni}$, which can be taken as representation of their neighbors, show some early variance. However, by $t_{\mathrm{f}}$, and even by $t_{\mathrm{expl}}$, these differences are only a few percent. There is a consistent and stark discrepancy between the in situ and post-processed total mass of $\mathrm{{}^{44}Ti}$ (yellow lines) and, to a lesser extent, $\mathrm{{}^{4}He}$, found in all four models. In B12-WH07, for example, the value from the in situ calculation, $M_{\mathrm{C}}^{\mathrm{\mathrm{{}^{44}Ti}}}(t_{\mathrm{f}})\approx 1.08\times 10^{-3}~{}M_{\sun}$, is greater than that from post-processing, $M_{\mathrm{PP}}^{\mathrm{\mathrm{{}^{44}Ti}}}(t_{\mathrm{f}})\approx 1.24\times 10^{-4}~{}M_{\sun}$, by nearly an order of magnitude. To understand the origin of these inconsistencies, consider that the nuclei most affected are products of $\alpha$-rich freeze-out occurring in low-density, expanding ejecta and are, therefore, most susceptible to inadequate tracer particle spatial resolution. The deficit of $\alpha$-rich freeze-out products in the post-processing results cannot be attributed to a lack of mixing therein, as mixing would dilute the $\alpha$-richness of this ejecta and, therefore, further reduce the production of these nuclei.

To provide context regarding the magnitude of this disagreement, we compare $M_{\mathrm{C}}^{i}(t_{\mathrm{f}})$ and $M_{\mathrm{PP}}^{i}(t_{\mathrm{f}})$ to post-processing calculations performed with a more realistic 150-species network in Figure 13. In each model, the unbound $\mathrm{{}^{56}Ni}$ mass is relatively unaffected, but there are significant differences in the total unbound mass for $\mathrm{{}^{44}Ti}$, $\mathrm{{}^{48}Cr}$, $\mathrm{{}^{52}Fe}$, and $\mathrm{{}^{60}Zn}$. These differences can be largely attributed to the availability of additional reaction pathways during explosive burning, particularly those involving $(n,\gamma)$ and $(p,\gamma)$ reactions (Hix & Thielemann 1996; Timmes et al. 2000; Magkotsios et al. 2010). The effect of spatial resolution on the production of $\mathrm{{}^{44}Ti}$ is of roughly the same order, but in the opposing direction, as that seen when switching from the $\alpha$-network to sn150. Without an in situ large network simulation, we are unable to fully quantify how spatial resolution of the tracer particles may impact nucleosynthesis with realistic nuclear networks. However, this work suggests that much higher tracer particle spatial resolution, and ultimately, large, in situ nuclear networks are needed to correctly calculate $\alpha$-rich freeze-out.

The susceptibility of post-processing to inadequate spatial resolution can also be broadly characterized by the mass distribution of $Y_{\mathrm{e}}(t_{\mathrm{f}})$ for unbound matter (see Figure 14). In these mass histograms, calculated from both individual zone data (blue) and tracer particle data (red), it becomes clear that the masses represented by a single tracer particle (shaded region) in the B-Series models fail to adequately resolve ejecta outside of $0.49<Y_{\mathrm{e}}(t_{\mathrm{f}})<0.51$. While this constitutes by far the majority of the matter, it neglects some of the most interesting nucleosynthesis. For B12-WH07 and, to a large extent, B15-WH07, increasing the number of tracer particles by an order of magnitude would serve to capture much more of the $Y_{\mathrm{e}}$ mass distribution and is a relatively easy solution to implement. Applying a similar approach to B20-WH07 and B25-WH07 would require an untenable number of particles to effectively capture the long tails of the distribution.

As in Equation 12, we define

[TABLE]

From this measure, it becomes clear that uncertainties stemming from spatial resolution, even though they apply to only a small subset of species, are the largest, single source of nucleosynthesis error in the B-Series models (see Table 2).

4.4 NSE transition and network size

Compositional evolution via the nuclear reaction network, both within the simulation and as a post-processing calculation, is an initial-value problem, with the initial values provided by one of two different methods (Section 3.4). For particles which have never reached NSE, the initial composition of the progenitor is mapped onto the nuclear network used for the nuclear evolution, providing the initial abundances. For particles which have reached NSE, the initial abundances from a NSE calculation can be mapped onto the network. Chimera’s treatment of the transition of matter into NSE (Section 3.2) is comparable to (or in some cases, better than) that used in other CCSN codes of similar capability—e.g., Castro (Almgren et al. 2010; Zhang et al. 2011, 2013), Prometheus-Vertex (Buras et al. 2006b; Marek & Janka 2009), CoCoNuT-Vertex (Müller et al. 2012b), Zelmani (Ott et al. 2011, 2013), and Zeus+IDSA (Suwa et al. 2010, 2013). The transition condition is motivated by the temperatures and densities at which complete silicon burning would occur within the current global timestep. For temperatures above this threshold, the use of the nuclear network is superfluous, as the network will achieve NSE every timestep. For simplicity, the transition out of NSE occurs when the temperature drops below this NSE condition (Equation 1 for Chimera). However, for the rapidly changing conditions in expanding CCSN matter, the assumption of NSE has been shown to break down when the temperature falls below 6 GK (Meyer et al. 1998).

This leaves a dilemma for this or any similar post-processing study. Is it better to be consistent with the NSE-to-network transition used within the supernova simulation, or should an earlier (higher temperature) transition out of NSE be adopted for the network? As a test of the NSE transition criteria in the B-Series simulations, we post-process the nucleosynthesis using the same $\alpha$-network used with the simulations and vary the conditions at which the transition to nuclear burning from the NSE composition occurs. In Figure 15, we compare the unbound masses of individual isotopes, transitioned out of NSE using either $T_{\mathrm{NSE}}(\rho)$, as defined in Equation 1, or $T_{\mathrm{NSE}}=8$ GK and show the relative deviation for each species $i$ as $\delta_{\mathrm{NSE}}^{i}\equiv\log_{10}(M_{\mathrm{PP}}^{i}(T_{\mathrm{NSE}}(\rho))/M_{\mathrm{PP}}^{i}(T_{\mathrm{NSE}}=8~{}\mathrm{GK}))$. The in situ calculation transitions out of NSE at a temperature $\approx$2–3 GK lower than the tested value of 8 GK. This reduces the $\alpha$-richness of the eventual freeze-out by maintaining NSE longer than it would physically. As a result, we see a significant shift in the masses of $\mathrm{{}^{4}He}$, $\mathrm{{}^{44}Ti}$, $\mathrm{{}^{48}Cr}$, and $\mathrm{{}^{60}Zn}$ in each model (e.g., $M_{\mathrm{PP}}^{\mathrm{\mathrm{{}^{44}Ti}}}(T_{\mathrm{NSE}}=8~{}\mathrm{GK})/M_{\mathrm{PP}}^{\mathrm{\mathrm{{}^{44}Ti}}}(T_{\mathrm{NSE}}(\rho))\approx 3.0$), highlighting the failure of $T_{\mathrm{NSE}}(\rho)$ to fully capture $\alpha$-rich freeze-out in the ejecta. This argues that for a stricter set of criteria for the breakdown of NSE should be adopted for use in post-processing than the $\approx$5–6 GK commonly used within models of the CCSN mechanism so that one may properly achieve $\alpha$-rich freeze-out. Furthermore, the transition used within the models themselves should be questioned if direct nucleosynthesis results are to be used to constrain models. We define

[TABLE]

for direct comparison with the uncertainties defined by Equation 12 and Equation 13 (see Table 2).

The same post-processing calculations using Chimera’s $\alpha$-network also provide insight into the limitations of this network on the nucleosynthesis in the simulation. Limitations in tracer particle spatial resolution significantly constrain our ability to capture $\alpha$-rich freeze-out in post-processing calculations (Section 4.3), but this uncertainty is mildly abated by improved treatment of the NSE transition, which enhances production of $\alpha$-rich freeze-out species (see Figure 15). This effect is visible for both sn150 and the $\alpha$-network, and $\delta_{\mathrm{NSE}}^{i}(\alpha)$ is nearly identical to $\delta_{\mathrm{NSE}}^{i}(\texttt{sn150})$. Differences in ejected mass of individual isotopes incurred from nuclear network size can vary (compare blue diamonds to red diamonds in Figures 15); the extent of which can be as high as an order of magnitude for $\mathrm{{}^{44}Ti}$. However, for species like $\mathrm{{}^{56}Ni}$, the effect is smaller ($\approx$20%). With this uncertainty in mind, we have reasonable confidence quoting the production of Ni and species from O–Ca from the in situ $\alpha$-network, but the values for the products of $\alpha$-rich freeze-out require the use of a larger in situ network.

4.5 Time resolution

The convective and turbulent nature of CCSNe make it highly probable for tracer particle thermodynamic states to change very rapidly. To capture this detail, in the B-series models, Chimera recorded the temporal history of each particle independently whenever the temperature changed by 0.1% from the last record and, likewise, limited changes in density and integrated neutrino number flux, $\phi_{\mathrm{\nu}}$, to 1%. These sampling criteria ensure that all significant features in the temporal histories of the tracer particles are captured. However, this high temporal cadence limits the number of tracers that can be evolved if the total cost of the simulation is not to be limited by tracer I/O. Independent sampling of individual tracers for the hundreds-of-thousands of tracer particles necessary to adequately sample a 3D simulation poses an even heavier and perhaps untenable load on I/O. One alternative to this approach is to record the thermodynamic states of all tracer particles at fixed time intervals along with the typical checkpoint data of the simulation. While this addresses the logistical concerns of managing a large amount of data, it runs the risk of insufficiently sampling the histories of individual particles. The biggest danger of less frequent sampling is the underestimate of local maxima in temperature and density in the thermodynamic history. The exponential temperature dependence of thermonuclear reactions can potentially greatly magnify the impact of under-sampling in the local peak conditions. The balanced reactions which maintain NSE, particularly those linking $\mathrm{{}^{4}He}$ to $\mathrm{{}^{12}C}$, are especially susceptible to peak temperature and density conditions (Meyer et al. 1998). Thermodynamic histories which fail to capture this peak behavior due to inadequate time resolution can thereby misestimate the conditions of $\alpha$-rich freeze-out.

We quantify this effect and other sampling-related concerns in tracer particle histories by post-processing thermodynamic profiles down-sampled to fixed time intervals. We compare the resulting composition profile to that from the dynamic time interval criteria described above for $\mathrm{{}^{4}He}$, $\mathrm{{}^{28}Si}$, $\mathrm{{}^{44}Ti}$, and $\mathrm{{}^{56}Ni}$ for both the $\alpha$-network and sn150 in Figure 16 for B12-WH07.

We find that a fixed rate of sampling the thermodynamic state reproduces post-processed abundances from the dynamic time interval criteria to within 0.5% so long as the fixed interval, $\delta t$, does not exceed 1 ms. For $\delta t>1$ ms, post-processing calculations fail to capture the important fluctuations in density and temperature, particularly those that define the conditions of $\alpha$-rich freeze-out. The subsequent errors incurred in the total unbound mass is independent of the reaction network size, and $\mathrm{{}^{44}Ti}$ is most affected by poor time resolution ($\approx$10% error for $\delta t=10$ ms and $\approx$30% error for $\delta t=30$ ms).

5 Gauging the completeness of the thermonuclear evolution

The goal of any calculation of CCSN nucleosynthesis is the final distribution of the ejecta in composition and velocity. The evolution to this final state can be divided into four rough stages: (1) the development of the explosion; (2) the development of the mass-cut that separates ejecta from proto-NS; (3) the completion of the nucleosynthesis as all ejecta expands and cools sufficiently for nuclear reactions to cease; and (4) the development of the final velocity distribution for the outgoing ejecta as the shock passes through the stellar envelope. In practice, these four stages are not so clearly demarcated as this list implies, and the overlap of these stages is exacerbated in multidimensional simulations.

For investigations primarily focused on the development of the explosion, there is a strong temptation to conserve computational resources by stopping the simulation once the successful revival of the supernova shock seems guaranteed or, at the latest, as the growth of the explosion energy levels off. From the nucleosynthesis perspective, such an abbreviation of the simulations certainly prevents examination of the final stage in the prior list; however, it is unclear how the intermediate stages may be affected. In this section, we attempt to quantify the impact of early termination of the simulation in terms of the uncertainties discussed in Section 4.

From the global analysis of the nuclear composition discussed in Bruenn et al. (2016, see Figure 22), it is clear that the nuclear composition is changing dramatically even 600 ms after bounce in all four models we consider here. This corresponds to the epoch when the explosion energy in these models begins to level off. We define this time as $t_{\mathrm{expl}}$ and estimate $t_{\mathrm{expl}}=600$ ms for B12-WH07 and $t_{\mathrm{expl}}=800$ ms for the other B-Series simulations (see Figure 12a in Bruenn et al. 2016). The decline in the growth of the explosion energy is concurrent with a similar decrease in the growth of the unbound mass ($M_{\mathrm{unb}}$; see Figure 1); however, using $\epsilon(t)$ in these figures as a proxy for the evolution of the multidimensional mass-cut, we see that the models at $t_{\mathrm{expl}}$ are much less mature than at $t_{\mathrm{f}}$. The disparity between $\epsilon(t_{\mathrm{expl}})$ and $[{\epsilon}](t_{\mathrm{expl}})$ clearly demonstrates the particles representative of ejecta are in flux, further supported by the growth in $M_{\mathrm{bound}}$ seen after $t_{\mathrm{expl}}$ in each model as the shock moves outward before leveling off near $t_{\mathrm{f}}$.

Such early truncation of the simulation often occurs with much of the ejecta at final temperatures above 3 GK (see Figure 6; yellow) and fails to capture all stages of explosive nucleosynthesis in situ when coupling to the hydrodynamics is important. In this case, investigation of the nucleosynthesis will much more heavily rely on extrapolations of the thermodynamic history, described in detail in Section 3.3.1. Coming, as it does, at a point the supernova’s evolution where hydrodynamics remains very active, this extrapolation is fraught with uncertainty. As a measure of this, we compare predictions of the nucleosynthetic yields determined from $t_{\mathrm{expl}}$ and $t_{\mathrm{f}}$ at 100 s in Figure 17. Whereas uncertainties relating to thermodynamic extrapolation at $t_{\mathrm{f}}$ (see Section 4.2) are small relative to other uncertainties, the magnitude of variations between $M_{+}^{i}(t_{\mathrm{expl}})$ and $M_{+}^{i}(t_{\mathrm{f}})$ for $A\geq 12$, which includes the effects of both a different multidimensional mass-cut and contrasting expansion timescales, can exceed an order of magnitude, with a factor of three being quite common. Interestingly, with the exception of $\mathrm{{}^{44}Ti}$, the $\alpha$-nuclei from silicon to nickel are relatively unaffected (e.g. $M_{+}^{\mathrm{\mathrm{{}^{56}Ni}}}(t_{\mathrm{expl}})/M_{+}^{\mathrm{\mathrm{{}^{56}Ni}}}(t_{\mathrm{f}})=1.033$). The production of these nuclei is largely tied to complete, explosive silicon burning ignited by the passage of the shock (see Table III in Woosley et al. 2002). As seen in Figure 6, by $t_{\mathrm{expl}}$, the shock is no longer heating newly swept up matter sufficiently ($T\gtrsim 5$ GK).

6 Summary

In this paper, we have discussed some of the uncertainties which complicate post-processing nucleosynthesis calculations from ab initio multidimensional CCSN simulations evolved beyond the initial stages of explosive nuclear burning but not yet to conditions for nuclear reactions to cease. We provide specific examples of how these uncertainties impact nucleosynthesis predictions for the four axisymmetric models of Bruenn et al. (2016). A detailed exposition of nucleosynthesis for these models using very large networks constitutes the subject of a forthcoming publication (Harris et al. 2017, in prep).

Our results can be summarized as follows:

Even after 1.2–1.4 s of post-bounce evolution and with asymptotic explosion energies, the multidimensional mass-cut remains unresolved in each model, impacting the production of nuclear species in borderline ejecta near ongoing accretion flows. The consequences of this are particularly pronounced in B25-WH07, wherein the ultimate fate of $\approx$0.2 $M_{\sun}$ of a total $\approx$1.3 $M_{\sun}$ of gravitationally-unbound matter is indeterminate. The state of the multidimensional mass-cut is even more dire if the simulations are truncated at 600–800 ms after bounce, after the explosion energy has begun to level off but before explosive nucleosynthesis has completed. 2. 2.

Despite temperatures below the threshold for explosive nuclear burning, a result of the extended running times, local hydrodynamic deviations from isentropic expansion continue to play a non-trivial role on secondary nuclear processes and neutrino-induced nucleosynthesis by altering the expansion timescale estimate necessary for extrapolation to freeze-out at $T\approx 0.5$ GK. 3. 3.

The ability of Lagrangian tracer particles to effectively reproduce the in situ nuclear burning conditions of the B-Series simulations is significantly reduced in regions of $\alpha$-rich freeze-out. In each of the four models, $\mathrm{{}^{44}Ti}$ is consistently under-produced. The magnitude of this effect is similar, but in the opposing direction, to that of replacing the $\alpha$-network used in the post-processing calculation with a more realistic 150-species network. While this effect targets only those species that result from $\alpha$-rich freeze-out, it has the largest global effect of the uncertainties tested. Ideally, large in situ networks are the best method to address this issue. Failing that, much higher numbers of tracer particles are needed. 4. 4.

Furthermore, we argue for a stricter set of criteria for transitioning out of NSE than the commonly used $T_{\mathrm{NSE}}\approx 5\textrm{--}6$ GK within models of the CCSN mechanism. While sufficient for the transition of matter into NSE, this criteria fails to capture the process of $\alpha$-rich freeze-out crucial to the production of $\mathrm{{}^{44}Ti}$. 5. 5.

We find that recording the thermodynamic conditions of all tracer particles at a fixed time interval $\delta t\leq 1$ ms to be a viable alternative to independently tracking the detailed history of each particle, limited by deviations in the nuclear burning inputs (i.e. $\rho$, $T$, and $\phi_{\mathrm{\nu}}$).

In light of these findings, we have modified Chimera, when possible, in an attempt to reduce these uncertainties for future models. Regarding Item 3, ongoing models improve the spatial resolution of tracer particle sampling by increasing the number of tracer particles by roughly one order of magnitude. The larger burden this would place on the file system is sufficiently alleviated with fixed time interval I/O for tracer particles. Our results in Section 4.5 (Item 5) give us confidence in this approach. Lastly, taking our own advice regarding Item 4 (Section 4.4), we have implemented a framework in Chimera to use arbitrary criteria for the NSE transition. When combined with quasi-statistical equilibrium methods (Hix et al. 2007; Parete-Koon et al. 2008), Chimera will be capable of seamlessly evolving nuclear burning networks through the NSE transition.

In theory, both the indeterminate mass-cut and expansion timescale uncertainties could be reduced by extending the simulation to freeze-out. However, given the inadequate spatial resolution of the tracer particles and an inherent limitation in the accuracy of the rate of nuclear energy released by the smaller network within the hydrodynamics, we cannot rely entirely on post-processing methods to obtain an accurate representation of the nucleosynthesis. Since the nucleosynthesis depends on the thermodynamic conditions and, consequently, nuclear energy generation, a feedback exists that cannot be captured with post-processing, significantly affecting the abundances of species such as $\mathrm{{}^{44}Ti}$, $\mathrm{{}^{57}Fe}$, $\mathrm{{}^{58}Ni}$, and $\mathrm{{}^{60}Zn}$ (Woosley & Weaver 1995). Improving upon the existing in situ $\alpha$-network with a more realistic 150-species nuclear network capable of properly tracking neutronization and energy release via particle captures is an important step towards resolving this problem and is the subject of ongoing work.

The research presented here was supported by the U.S. Department of Energy, Office of Science, Offices of Nuclear Physics and Advanced Scientific Computing Research, the NASA Astrophysics Theory and Fundamental Physics Program (NNH11AQ72I) and the National Science Foundation Nuclear Theory Program (PHY-1516197). The simulations described herein were performed via NSF TeraGrid resources provided by the National Institute for Computational Sciences under grant number TG-MCA08X010; resources of the National Energy Research Scientific Computing Center, supported by the U.S. DOE Office of Science under Contract No. DE-AC02-05CH11231; and an award of computer time from the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program at the Oak Ridge National Leadership Computing Facility, supported by the U.S. DOE Office of Science under Contract No. DE-AC05-00OR22725.