Insights into the Production of ⁴⁴Ti and Nickel Isotopes in Core-Collapse Supernovae

Abstract

We report nucleosynthetic results for both ⁴⁴Ti and nickel isotopes for eighteen three-dimensional (3D) core-collapse supernova (CCSN) simulations extended to $\sim$20 seconds after bounce. We find that many of our long-term models are able to achieve ⁴⁴Ti/⁵⁶Ni ratios similar to that observed in Cassiopeia A, and modern supernova models can synthesize up to $2\times 10^{-4}M_{\odot}$ of ⁴⁴Ti. Neutrino-driven winds and the fact that there can be simultaneous accretion and explosion in 3D models of core-collapse supernovae play central roles in its production. We conclude that the ⁴⁴Ti underproduction problem in previous CCSN models is no longer an issue. In addition, we discuss the production of both ⁵⁷Ni and stable nickel/iron ratios and compare our results to observations of SN1987A and the Crab.

1 Introduction

Supernova remnants offer unique laboratories for the study of core-collapse supernova (CCSN) explosions. Various aspects of the CCSN theory can be constrained by measuring the abundance and distribution of different nuclei in the remnants. Despite of the fact that CCSNe are able to produce many of the elements in the periodic table (Arcones & Thielemann, 2023), there are a few isotopes (such as the radioactive isotopes ⁴⁴Ti, ⁵⁶Ni, and ⁵⁷Ni) that are of particular interest both theoretically and observationally. These isotopes are created in the core of the supernova explosion, and as a result their formation is closely related to the details of the explosion mechanism itself. Two of them play important roles in the evolution of supernovae light curves because they radioactively heat the ejecta. Furthermore, with relatively short half-lives such isotopes are ideal targets for $\gamma-$ray astronomy (Kurfess et al., 1992; Iyudin et al., 1994, 1998; Diehl & Timmes, 1998; Renaud et al., 2006; Boggs et al., 2015; Weinberger et al., 2020). In addition, supernovae produce interesting stable nickel and iron isotopes, either directly or after decay. Hence, understanding the production and distribution of such isotopes in core-collapse explosions offers us a direct connection between theory and observation. Given this, it is important to make reliable theoretical predictions of supernova nucleosynthetic yields.

It is well known that multi-dimensional effects play an important role in CCSN explosions (for reviews, see Janka (2012); Burrows & Vartanyan (2021); Wang et al. (2022)), and such effects can significantly influence nucleosynthetic results. For example, the high ⁴⁴Ti masses (at least $10^{-4}M_{\odot}$) observed in SN1987A and Cassiopeia A (Cas A) are strongly inconsistent with spherically-symmetric predictions (The et al., 2006; Sukhbold et al., 2016; Curtis et al., 2019; Sieverding et al., 2023a). To correctly capture the role of multi-dimensional effects and to make reliable theoretical predictions of CCSN ejecta abundances, high-fidelity three-dimensional state-of-the-art CCSN simulations are required. However, due to computational resource limitations, previous 3D CCSN models (Wongwathanarat et al., 2017; Müller et al., 2017; Stockinger et al., 2020; Bollig et al., 2021; Sandoval et al., 2021; Sieverding et al., 2023a) that generated nucleosynthetic yields had at least one of the following unsatisfactory aspects: (a) the neutrino-transport scheme was approximate and the ejecta electron fractions ($Y_{e}$) were very uncertain; (b) the radiation transport calculation (or the simulation itself) was turned off long before reaching the asymptotic stage of explosion and the contribution of neutrino-driven winds was ignored; (c) the calculation of multi-dimensional hydro/radiation effects was switched to sphericized and/or parameterized models after only one or two seconds post-bounce when asymmetric long-lasting accretion is still strong; and (d) the size of the nucleosynthesis network was too small to reliably estimate the production of many isotopes. Most importantly, since previous work focused on a small number of models, the robustness and the progenitor dependencies of the results remained unknown.

In this work, we report the nucleosynthetic results of 18 3D CCSN simulations using the code Fornax (Skinner et al., 2019). The simulations were evolved to late times after bounce (many seconds, see Table 1) with full neutrino transport and multi-dimensional effects, and the asymptotic explosion energies have been reached for most models. By extrapolating both the tracer thermal histories and the neutrino-driven wind mass outflow rates, we are able to calculate the ejecta isotopic abundances to about twenty seconds post-bounce. Hence, for the first time one is able to perform a systematic study of the asymptotic nucleosynthetic yields in 3D CCSN models, including the contribution of neutrino-driven winds.

The paper is arranged as follows: We summarize the simulation specifics and the approach to the nucleosynthetic calculations in Section 2. In this section we also describe the extrapolation methods adopted. The main results are presented in Section 3, in which we provide a detailed discussion of the production of the radioactive isotopes ⁴⁴Ti, ⁵⁶Ni, and ⁵⁷Ni. We find that all of our 3D models are able to reach a ⁴⁴Ti/⁵⁶Ni ratio similar to that of Cas A, and we find that many models can produce more than $10^{-4}M_{\odot}$ of ⁴⁴Ti. Therefore, the high ⁴⁴Ti problem (The et al., 2006) of the past is naturally solved by long-term self-consistent 3D CCSN models, without the need to introduce fast initial rotation or jet-like structures. We find that the ⁵⁷Ni/⁵⁶Ni ratio of SN1987A (though with large errorbars) is only about 1.5 times higher than our models. Since ⁵⁷Ni is a neutron-rich isotope, it is possible that SN1987A contains slightly more neutron-rich ejecta than expected, which might arise from a slightly more prompt explosion than our models would suggest. We study the influence of different extrapolation methods in Section 3.2 and we also briefly discuss the stable Ni/Fe ratios we find. Finally, in Section 4 we summarize the paper and discuss the uncertainties and caveats of our work.

2 Method

In this paper, we provide a subset of the nucleosynthetic results for eighteen CCSN models with solar-metallicity progenitors sporting ZAMS masses that range from 9 to 60 $M_{\odot}$. Two models (the 19.56 and 40 $M_{\odot}$) form black holes and experience successful explosions (Burrows et al., 2023, 2024). All these models have already been published in Wang & Burrows (2024a) and Burrows et al. (2024). Table 1 summarizes some basic properties of the models. Since the publication of Wang & Burrows (2024a) and Burrows et al. (2024), some of the models (the 17, 18, 18.5, 20, 25, 60 $M_{\odot}$ models) have been simulated to even later times after bounce and these updates are reflected in what we present here.

All models are generated using the multi-group radiation/hydrodynamics code Fornax (Skinner et al., 2019; Vartanyan et al., 2019a; Burrows et al., 2019, 2020) with 12 logarithmically-distributed energy groups for each of our three neutrino species (electron-type, anti-electron type, and “$\mu$-type,” a bundling of the $\mu$ and $\tau$ neutrinos and anti-neutrinos). The SFHo nuclear equation of state (EOS) of Steiner et al. (2013) is used. Each simulation is done with a spherical grid with (1024, 128, 256) cells along the $r$, $\theta$, and $\phi$ directions, and the outer boundary radii vary from 30,000 to 100,000 kilometers (km). Each model is first simulated spherically until 10 milliseconds after bounce, after which it is continued in full 3D. No perturbations are introduced when the 3D simulations are started, except in the 9a model. The initial perturbations in 9a model cause a slightly more prompt explosion with more neutron-rich ejecta, which is very different from the perturbation-free 9b model of the same progenitor (Wang & Burrows, 2024a). However, as shown in Wang & Burrows (2024a), realistic initial perturbations in more massive models don’t seem to influence the electron fraction ($Y_{e}$) or nucleosynthetic yields of the ejecta.

Nucleosynthetic yields in this paper are calculated in the same way as in Wang & Burrows (2024a). We use post-processed backwardly-integrated tracer particles to track the dynamical and thermal evolution of the ejected matter. About 320,000 tracers are injected into each simulation. The tracers are placed logarithmically along the r-direction above 500 km and uniformly along the $\theta$- and $\phi$-directions at the end of each simulation. A more detailed description of the tracer method can be found in Wang & Burrows (2024a). The tracer trajectories with neutrino field information are fed into SkyNet (Lippuner & Roberts, 2017), and a 1530-isotope network including elements up to $A=100$ is used to calculate the nucleosynthetic results. The reaction rates are taken from the JINA Reaclib (Cyburt et al., 2010) database, and we include neutrino interactions with protons and neutrons. However, reactions for the $\nu-$process (Woosley et al., 1990) are not included. The NSE (nuclear statistical equilibrium) criterion is set at 0.6 MeV ($\sim$7 GK). The electron fractions ($Y_{e}$) are calculated by Fornax when the temperature is above the NSE threshold, which allows the neutrino spectra to be appropriately non-thermal. A nucleosynthesis calculation starts from the point after which the tracers never reach NSE again. We use the Lodders (2021) values whenever solar abundances are referred to in this paper.

Although our 3D CCSN models have been evolved to many seconds post-bounce (see Table 1), we find that some of them have not been carried long enough to reach the asymptotic abundances of certain isotopes. The first reason for this is that not all the tracers representing the ejecta have cooled down and finished their nucleosynthetic evolution, especially those ejected in the neutrino-driven wind phase. The second reason is that even at the end of our longest simulations, the neutrino-driven wind is still actively producing certain isotopes at a non-negligible rate. Therefore, to get the final ejecta isotopic abundance, two kinds of extrapolations are required: (a) extrapolation of the tracer thermal histories to calculate the asymptotic abundances of the matter ejected before the end of our simulations and (b) extrapolation of the mass outflow rates of the neutrino-driven winds to estimate how much more matter will be ejected after the end of our 3D simulations.

As a result, we extrapolate the tracer thermal histories of the 3D simulations using a power-law method (Ning et al., 2007; Zha et al., 2024). The temperatures and densities at the end of our simulations are continued using

$$\begin{split}T(t)=&T(t_{\rm end})(1+(t-t_{\rm end})/\tau)^{-2/3},\\ \rho(t)=&\rho(t_{\rm end})(1+(t-t_{\rm end})/\tau)^{-2},\\ \end{split}$$ (1)

where $t_{\rm end}$ is the end time of each simulation, and $\tau$ is the expansion timescale fitted from the last 10 ms of the associated trajectories of the 3D simulations. There is also an exponential method used in the literature for trajectory extrapolations (Magkotsios, 2011). The impact on nucleosynthetic results caused by different trajectory extrapolation methods will be discussed in more detail in Section 3.2.

Figure 1 shows the mass outflow rates of the neutrino-driven winds measured at 500 km as a function of time for all our models, smoothed over 100 milliseconds. We fit an exponential relation to the last one second of each simulation. We choose the exponential form not only because it fits the data well, but also to make sure that the total amount of mass ejected in the extrapolated phase is always convergent. We assume that the future ejected matter will have the same asymptotic abundances as the average values of tracers ejected in the last second. The black-hole formers don’t need this type of extrapolation because their neutrino-driven winds are instantly turned off after black hole formation (which is at the end of the 3D simulation).

With the extrapolation methods described above, we calculate the nucleosynthetic evolution of all models to about twenty seconds post-bounce. This allows us to study the final yields of the important isotopes we discuss in this work.

$M_{\rm ZAMS}$ [$M_{\odot}$]	$t_{\rm end}$ [s]	$E_{\rm exp}$ [B]	$M_{\rm NSE}$ ($M_{\rm extra})$[$M_{\odot}$]	44Ti [$10^{-4}M_{\odot}$]	56Ni [$M_{\odot}$]	57Ni [$10^{-3}M_{\odot}$]	Stable Ni/Fe
9(a)	1.775	0.111	0.0128 ($4.88\times 10^{-4}$)	0.013	0.0019	0.061	0.502
9(b)	1.950	0.094	0.0115 ($3.35\times 10^{-4}$)	0.035	0.0063	0.065	0.055
9.25	3.532	0.124	0.0198 ($3.17\times 10^{-4}$)	0.074	0.0106	0.116	0.054
9.5	2.375	0.142	0.0262 ($2.24\times 10^{-3}$)	0.121	0.0156	0.180	0.054
11	4.492	0.321	0.0618 ($3.02\times 10^{-2}$)	0.413	0.0418	0.763	0.216
15.01	4.383	0.288	0.0703 ($5.64\times 10^{-3}$)	0.680	0.0548	1.212	0.056
16	4.184	0.390	0.0823 ($1.24\times 10^{-2}$)	0.648	0.0653	1.025	0.045
17	6.164	1.123	0.1889 ($7.63\times 10^{-3}$)	1.443	0.1123	2.706	0.079
18	8.617	0.516	0.1480 ($2.44\times 10^{-2}$)	1.055	0.1181	2.219	0.054
18.5	6.408	1.155	0.1760 ($8.96\times 10^{-3}$)	1.192	0.1421	2.335	0.043
19	4.075	0.466	0.1118 ($7.06\times 10^{-2}$)	1.095	0.0949	1.515	0.052
19.56∗	3.890	2.451	0.3203 (0)	1.619	0.2536	3.594	0.047
20	6.429	0.881	0.1437 ($1.28\times 10^{-2}$)	1.080	0.1087	2.126	0.034
23	6.228	0.513	0.1223 ($1.92\times 10^{-2}$)	1.180	0.0944	1.627	0.055
24	3.919	0.779	0.1586 ($2.64\times 10^{-2}$)	1.287	0.1348	2.155	0.039
25	6.419	1.101	0.2419 ($3.18\times 10^{-2}$)	2.124	0.1887	3.128	0.037
40∗	1.760	1.804	0.2151 (0)	1.162	0.1798	4.661	0.106
60	8.026	0.541	0.1616 ($4.82\times 10^{-2}$)	1.645	0.1321	2.785	0.117

Table 1: This table summarizes some basic properties of our models. The core bounce time is set to $t=0$. $t_{\rm end}$ is the time after bounce when the full 3D simulation stops, $E_{\rm exp}$ is the explosion energy in Bethes ($\equiv$10⁵¹ ergs), $M_{\rm NSE}$ is the mass of ejecta that has ever reached our NSE threshold (0.6 MeV), and $M_{\rm extra}$ is the mass of the extrapolated neutrino-driven winds. The isotopic yields are calculated at about twenty seconds post-bounce. The solar stable Ni/Fe ratio is 0.058 (Lodders, 2021), and the ratios in our models (including the residual progenitor envelope) range between about 0.6 to 7 times the solar value. The 19.56 and 40 $M_{\odot}$ models form black holes at the end of the simulations (Burrows et al., 2023, 2024); thus, they have $M_{\rm extra}=0$.

Figure 1: Neutrino-driven wind mass outflow rates measured at 500 km, smoothed over 100 milliseconds. Although there are fluctuations in the wind mass outflow rates, they generally decay exponentially with time at later times, which allows us to extrapolate their evolution. We fit exponential relations to the wind mass outflow rates in the last 1.5$-$2 seconds of each simulation. The black-hole formers (19.56 and 40 $M_{\odot}$) are not shown on this plot because their neutrino-driven winds are instantly turned off after the black hole formation (which is the end of the 3D simulation).

Figure 2: ⁵⁶Ni yield as a function of time. The circular dot on each curve marks the beginning point of extrapolation. The black-hole formers (19.56 and 40 $M_{\odot}$) are plotted as dashed curves, and they have no extrapolated neutrino-driven winds. In all models, the majority of ⁵⁶Ni is produced in the first 4$-$5 seconds after bounce, while the more aspherical models with stronger long-lasting accretion (e.g., the 11, 19, 25, and 60 $M_{\odot}$ models) experiences a further 10% ⁵⁶Ni production during the extrapolated phase. This late-time ⁵⁶Ni increase results from the $\alpha$-rich freeze-out process in the enhanced neutrino-driven winds (see Figure 1) due to their longer-lasting accretion. However, in general after a few seconds the neutrino-driven winds don’t contribute much to the total ⁵⁶Ni yields.

Figure 3: ⁴⁴Ti yield (left) and ⁴⁴Ti/⁵⁶Ni ratio (right) and as a function of time. The circular dot on each curve marks the beginning point of extrapolation. The black-hole formers (19.56 and 40 $M_{\odot}$) are plotted as dashed curves. The black horizontal line shows the solar ⁴⁴Ca/⁵⁶Fe ratio as a reference. The asymptotic ⁴⁴Ti yield is not achieved in most models until 10 seconds post bounce, except the lowest-mass ones (e.g. 9a, 9b, 9.25, and 9.5 $M_{\odot}$). The growth of the ⁴⁴Ti yield in the extrapolated phase comes mostly from the $\alpha$-rich freeze-out process in the proton-rich neutrino-driven winds, which initially produces proton-rich isotopes such as ⁴⁵V and ⁴⁶Cr and then converts them into ⁴⁴Ti as the temperature decreases. See text for a discussion.

Figure 4: Final ⁴⁴Ti and ⁵⁶Ni yields of our theoretical 3D models compared to observed values in Cas A and SN1987A. The black-hole formers (19.56 and 40 $M_{\odot}$) are plotted as square dots. In addition to our 3D models, we include the 18.88 $M_{\odot}$ model from Sieverding et al. (2023a) (the black triangle, labeled as S23). The observed ⁴⁴Ti and ⁵⁶Ni masses of Cas A are taken from Grefenstette et al. (2017) and Hwang & Laming (2012), while the SN1987A values are taken from Boggs et al. (2015) and McCray & Fransson (2016). Our 3D models can easily produce the high amount of ⁴⁴Ti and reach the high ⁴⁴Ti/⁵⁶Ni ratio observed in Cas A. Therefore, the “⁴⁴Ti-problem” (The et al., 2006) is no longer an issue in self-consistent long-term, 3D, initially non-rotating CCSN models. In addition, the lower ⁴⁴Ti/⁵⁶Ni ratios of the black-hole formers clearly show the important role of neutrino-driven winds in the production of ⁴⁴Ti, since the winds in these models are turned off after black hole formation.

Figure 5: ⁵⁷Ni yield and ⁵⁷Ni/⁵⁶Ni ratio as a function of time. The circular dot on each curve marks the beginning point of extrapolation. The black-hole formers (19.56 and 40 $M_{\odot}$) are plotted as dashed curves. The black horizontal line shows the solar ⁵⁷Fe/⁵⁶Fe ratio as a reference. Similar to ⁴⁴Ti, a large amount of ⁵⁷Ni is produced at later times ($>$5 seconds post-bounce) due to the $\alpha$-rich freeze-out process in the proton-rich neutrino-driven winds. However, a high ⁵⁷Ni/⁵⁶Ni ratio (e.g., $3\times 10^{-2}$) can’t be reached unless the explosion has a significant amount of neutron-rich ejecta.

Figure 6: Final ⁵⁷Ni and ⁵⁶Ni yield of 3D theoretical models compared to observed values in SN1987A. The black hole formers (19.56 and 40) are plotted as square dots. Despite of the large error bars, SN1987A has a higher ⁵⁷Ni/⁵⁶Ni ratio than most of our models. Only the most neutron-rich explosion (9a) shows a ⁵⁷Ni/⁵⁶Ni ratio within the one-sigma range of the observation. One possibility is that SN1987A exploded more prompt than most of our massive models and ejected a larger amount of neutron-rich matter. See text for more discussions.

3 Results

3.1 ⁴⁴Ti, ⁵⁶Ni, and ⁵⁷Ni

⁴⁴Ti, ⁵⁶Ni, and ⁵⁷Ni are produced in CCSN explosions in measurable amounts and the radioactive decay of the latter two can help power supernova light curves. In addition, the presence and yields of these isotopes can be inferred by direct $\gamma$-ray and X-ray observations.

A range of observations have indicated that in both SN1987A and Cassiopeia A (Cas A) at least $10^{-4}M_{\odot}$ of ⁴⁴Ti was produced, with ⁴⁴Ti/⁵⁶Ni ratios above $10^{-3}$ (Hwang et al., 2004; Hwang & Laming, 2012; Grefenstette et al., 2014; Boggs et al., 2015; McCray & Fransson, 2016; Grefenstette et al., 2017). These values were inconsistent with previous theoretical CCSN models, in particular those that were spherically-symmetric (The et al., 2006; Sukhbold et al., 2016; Curtis et al., 2019; Sieverding et al., 2023a) and this ⁴⁴Ti problem was identified almost two decades ago (The et al., 2006). One possible solution is the effect of explosion asymmetry. The strongly aspherical, though axisymmetric (2D), models of Nagataki et al. (1997) and Maeda & Nomoto (2003) were able to obtain high ratios of ⁴⁴Ti/⁵⁶Ni. However, these models used parameterized initial conditions and didn’t explain the origin of the strong, even jet-like, asymmetries necessary. Such jet-like explosions are not commonly seen in self-consistent 3D simulations (Lentz et al., 2015; Wongwathanarat et al., 2015; Roberts et al., 2016; Wongwathanarat et al., 2017; Müller et al., 2017; O’Connor & Couch, 2018; Vartanyan et al., 2019b; Ott et al., 2018; Summa et al., 2018; Glas et al., 2019; Burrows et al., 2019; Nagakura et al., 2019; Burrows et al., 2020; Stockinger et al., 2020; Bollig et al., 2021; Sandoval et al., 2021; Burrows et al., 2024), unless the progenitor is initially rotating fast, which is thought to be a rare occurrence (Heger et al., 2005; Obergaulinger & Aloy, 2021; Powell et al., 2023; Shibagaki et al., 2024). Three-dimensional simulations done by Wongwathanarat et al. (2017) could explain the ⁴⁴Ti production in Cas A, but those authors used prescribed neutrino luminosities at the inner grid boundary to tune the explosion energy to a desired value. The gray neutrino transport treatment they used was approximate, and their models could not produce enough ⁴⁴Ti if the tracer electron fractions ($Y_{e}$s) were taken from the 3D simulations instead of by manually setting them to a value of 0.5.

Recently, Sieverding et al. (2023a) published a self-consistent 3D initially non-rotating CCSN model that reached a ⁴⁴Ti/⁵⁶Ni ratio close to observations. This paper emphasized the importance of long-term 3D simulations because (a) tracer trajectories in 3D show more complex temperature and density evolution and (b) the total ejecta mass continues to grow for many seconds. Both effects can aid in the production of ⁴⁴Ti. Although this work is an important step in solving the ⁴⁴Ti problem, the total amount of ⁴⁴Ti in their model is still below $10^{-4}M_{\odot}$ and does not explain the observed values in Cas A and SN1987A. It was also unclear from their paper if the high ⁴⁴Ti/⁵⁶Ni ratio is progenitor model dependent. Therefore, a systematic study of CCSN ⁴⁴Ti production with many long-term 3D simulations has been lacking. In this section, we present for the first time such a systematic study of ⁴⁴Ti production in initially non-rotating 3D CCSN explosion models, together with a corresponding discussion of ⁵⁶Ni and ⁵⁷Ni production and the effects of neutrino-driven winds.

We start with a discussion of the production of ⁵⁶Ni. In CCSN explosions, ⁵⁶Ni is produced explosively by the shock wave as it traverses the mantle matter (predominantly the inner reaches of the oxygen shell) and in the neutrino-driven winds via $\alpha$-rich freeze-out (Woosley et al., 2002; Nomoto et al., 2013; Wanajo, 2023; Arcones & Thielemann, 2023). Explosive nucleosynthesis terminates quickly as the shock expands and the post-shock temperature drops, while the neutrino-driven winds can last many seconds post-bounce. Wang & Burrows (2023) have shown that long-lasting accretion can strengthen the neutrino-driven winds by enhancing the neutrino luminosities and increasing the density in the wind formation region. Therefore, neutrino-driven winds play a more important role in nucleosynthesis in more aspherical 3D explosions with stronger long-lasting accretion (usually the case in more massive models) than in 1D spherical models.

Figure 2 shows the temporal evolution of ⁵⁶Ni production in each simulation. The circular dot marks the point when model extrapolation begins. In all models, the vast bulk of ⁵⁶Ni is produced in the first 4$-$5 seconds after bounce, meaning that the explosive nucleosynthesis component and the early-phase neutrino-driven winds dominate ⁵⁶Ni production. More aspherical models with stronger long-lasting accretion (e.g., the 11, 19, 25, and 60 $M_{\odot}$ models) experience only $10-20$% extra ⁵⁶Ni production during the extrapolated phase. This late-time ⁵⁶Ni growth comes from the $\alpha$-rich freeze-out process in the enhanced neutrino-driven winds due to long-lasting accretion (see Figure 1). The wind mass outflow rates of these models have shallower slopes and, thus, last for a longer period of time. However, in general after a few seconds the neutrino-driven winds don’t contribute much to the total ⁵⁶Ni yields.

The production of ⁴⁴Ti is different from that of ⁵⁶Ni. Magkotsios (2011) has studied in detail the production channels of ⁴⁴Ti, and we have used this work to identify the most important channel for ⁴⁴Ti production. We find that the dominant contribution to ⁴⁴Ti production in our simulations comes from the neutrino-driven winds and corresponds to the “(p,$\gamma$)-leakage” regime discussed in Magkotsios (2011). This is a special type of freeze-out. Briefly, in a proton-rich environment (which is very common in neutrino-driven winds), the (p,$\gamma$) reactions transfer the nucleosynthetic flow from symmetric to proton-rich isotopes. As the temperature decreases, the freeze-out process starts and the weak interactions decrease $Y_{e}$ and cause the flow to transfer back towards more stable isotopes. In the case of ⁴⁴Ti, this means that isotopes such as ⁴⁵V and ⁴⁶Cr are produced at early times, and they are converted into ⁴⁴Ti as the temperature decreases. This explains the non-monotonic growth of the ⁴⁴Ti abundances in our tracers. The decay of ⁴⁴V to ⁴⁴Ti also contributes to ⁴⁴Ti growth, but plays only a secondary role. The freeze-out process remains active even when the temperature is lower than 1 GK. Therefore, the long production timescale of ⁴⁴Ti (Sieverding et al., 2023a; Wang & Burrows, 2024a) is due to two main factors: (a) the neutrino-driven winds last for many seconds and keep contributing to the ⁴⁴Ti yield and (b) the freeze-out process requires a very low temperature to reach the asymptotic state, and this takes a long time to achieve during the expansion. Figure 3 shows the ⁴⁴Ti yields and ⁴⁴Ti/⁵⁶Ni ratios of all our models as a function of time. The circular dot on each curve marks the point when the extrapolation is started. Many models take more than 10 seconds to reach their asymptotic ⁴⁴Ti yields, and $10^{-4}M_{\odot}$ of ⁴⁴Ti production can easily be achieved by models more massive than $\sim$16 $M_{\odot}$. More massive models in general produce more ⁴⁴Ti (see Table 1). The right panel shows that the more massive models also tend to have higher ⁴⁴Ti/⁵⁶Ni ratios than the low mass models (e.g., 9a, 9b, 9.25, and 9.5 $M_{\odot}$), so the more massive models not only have higher ejecta masses but also a higher ⁴⁴Ti production efficiency. This is because more massive explosion models are in general more asymmetric than the low-mass progenitor models due to their generally more massive envelopes.

Multi-dimensional effects play a central role in the production of ⁴⁴Ti. As Sieverding et al. (2023a) has pointed out, large-scale convection leads to non-monotonic thermal histories. Thus, nucleosynthesis can experience multiple nucleosynthetic phases, which aid the production of ⁴⁴Ti. We see the same behavior in our 3D simulations. However, this effect can’t solely explain the significantly higher ⁴⁴Ti yields in our simulations, since the ejecta with significantly non-monotonic thermal histories contribute only a small fraction to the ⁴⁴Ti yields. The most important 3D effect is the so-called “simultaneous explosion and accretion” phenomenon. This phenomenon is sometimes referred to as “long-lasting” or “fallback” accretion. In the context of this phenomenon, the stellar envelope above the proto-neutron star (PNS) is not immediately swept away by the shock wave after the launch of the explosion. Instead, it forms a few inward-directed funnels that extend almost to the PNS surface. Matter in the outer layers can fall very close to the PNS through these funnels, after which it is heated by neutrinos and ejected as part of the anisotropic neutrino-driven wind. The simultaneous explosion and accretion phenomenon significantly increases the amount of mass interacting with neutrinos, which adds to the $\alpha$-rich freeze-out component. This enhancement of $\alpha$-rich freeze-out lasts for many seconds, until accretion is completely terminated and a spherical wind develops. In addition, interactions with neutrinos push $Y_{e}$ to higher values, which provides the proton-rich environment required by the “(p,$\gamma$)-leakage” channel. Therefore, the simultaneous explosion and accretion phenomenon (redirecting fallback into a component of the late-time winds) aids the production of ⁴⁴Ti by enhancing the proton-rich neutrino-driven winds. This is the origin of the high ⁴⁴Ti yields we observed in our long-term 3D CCSN simulations.

Neutrino-driven winds are important for the production of ⁴⁴Ti because they are the major context of $\alpha$-rich freeze-out. Thus, it is crucial to capture the behavior of neutrino-driven winds with tracers. For post-processed tracers, the backward integration method has two major advantages: (a) this method inserts the tracers at the end of the simulations, so all ejecta regions can be identified and represented and (b) this method reproduces the freeze-out conditions more accurately than the forward method, since it does not need to follow the tracers through the complicated fluid motions in the neutrino-driven convection region (Reichert et al., 2023; Sieverding et al., 2023b; Wang & Burrows, 2023). Therefore, the backward method gives more accurate nucleosynthetic results than the forward method.

Figure 4 plots the ⁴⁴Ti and ⁵⁶Ni yields predicted by our 3D theoretical CCSN models, along with the observed Cas A and SN1987A yields. The observed ⁴⁴Ti and ⁵⁶Ni masses of Cas A are taken from Grefenstette et al. (2017) and Hwang & Laming (2012), while the SN1987A values are taken from Boggs et al. (2015) and McCray & Fransson (2016). Our 3D models can easily produce high amounts of ⁴⁴Ti and reach the high ⁴⁴Ti/⁵⁶Ni ratio observed in Cas A. In addition, the lower ⁴⁴Ti/⁵⁶Ni ratios of the black-hole formers clearly show the important role of neutrino-driven winds in the production of ⁴⁴Ti, since the winds in these models are turned off after the black hole formation. Therefore, the “⁴⁴Ti-problem” (The et al., 2006) is no longer an issue in self-consistent long-term 3D initially non-rotating CCSN models. SN1987A is still about two sigmas above the ⁴⁴Ti/⁵⁶Ni ratios predicted by our models. We now turn to a discussion of this issue and its corresponding ⁵⁷Ni/⁵⁶Ni ratio.

⁵⁷Ni is an interesting radioactive isotope. Figure 5 shows the temporal evolution of the ⁵⁷Ni yields and ⁵⁷Ni/⁵⁶Ni ratios in our models. Two classes of behaviors are seen. If the explosion has ever experienced a neutron-rich phase (e.g., models 9a, 11, and 17), the ⁵⁷Ni/⁵⁶Ni ratio quickly rises to a high value during the neutron-rich phase and stays there. Although ⁵⁶Ni production in the neutrino-driven winds may decrease this ratio a bit, the overall ⁵⁷Ni/⁵⁶Ni ratios of these models remain relatively high. For models with no neutron-rich ejecta, the production timescale of ⁵⁷Ni is even longer than that of ⁴⁴Ti. This is because ⁵⁷Ni is also produced in significant amounts in the $\alpha$-rich freeze-out process in proton-rich neutrino-driven winds. The asymptotic ⁵⁷Ni/⁵⁶Ni ratios in these two different classes are not as different as they are initially, unless the neutron-rich phase occurs so early that the explosive nucleosynthesis component (which has more mass than the winds) partly becomes neutron-rich (witness 9a).

Figure 6 shows the ⁵⁷Ni and ⁵⁶Ni yields predicted by 3D theoretical CCSN models and observed in SN1987A. The ⁵⁷Ni amount of SN1987A is taken from Kurfess et al. (1992) and the ⁵⁶Ni amount is taken from Hwang & Laming (2012). Similar to the ⁴⁴Ti case, SN1987A has a ⁵⁷Ni/⁵⁶Ni ratio about 1.5 sigmas above the values seen in our 3D simulations. One possible explanation for the higher ⁴⁴Ti/⁵⁶Ni and ⁵⁷Ni/⁵⁶Ni ratios in SN1987A is that the explosion might have had an early and strong neutron-rich phase, so that a fraction of the matter participating in explosive nucleosynthesis becomes neutron-rich. This will enhance ⁵⁷Ni production, but the most significant influence is that the production of ⁵⁶Ni at the early phase will be suppressed. ⁴⁴Ti is not influenced by this early neutron-rich phase since it is produced mostly in the wind phase. Therefore, a more neutron-rich explosion, which might be caused by a more prompt explosion, will naturally lead to higher ⁴⁴Ti/⁵⁶Ni and ⁵⁷Ni/⁵⁶Ni ratios. This is one possibility for SN1987A.

As the name indicates, ⁴He can be produced in significant amounts by the $\alpha$-rich freeze-out. Therefore, ⁴He exists both in the helium envelope and in the innermost ejecta where the ⁵⁶Ni abundance is high. The coexistence of ⁴He with ⁵⁶Ni is interesting because it can absorb the $\gamma$-rays emitted by the decay of ⁵⁶Ni and influence the emergent spectra. ⁴He in the nickel/iron-rich region of supernova remnants might be a potentially interesting observable. We find that the ⁴He/⁵⁶Ni ratios in the ⁵⁶Ni-rich ejecta ($X_{{}^{56}\rm Ni}>10^{-5}$) are around 0.4$-$0.5 in most our models, while the low-mass models (9a, 9b, 9.25, 9.5, and 11 $M_{\odot}$) have higher ⁴He/⁵⁶Ni ratios around 0.8 .

3.2 Effects of the Extrapolation Method

Since a significant amount of ⁴⁴Ti and ⁵⁷Ni is produced during the extrapolated phase, will the above results be changed by the choice of extrapolation method? The contribution from the extrapolated neutrino-driven wind accounts for no more than 20% of the yield of these radioactive isotopes, so the wind extrapolation methods should not influence the results to any significant degree. However, it is not impossible that how we chose to extrapolate the tracer thermal histories might change our conclusions, since this influences most of the ejecta that never reaches high temperatures. Therefore, it would be helpful to have a robustness test of our choice of tracer extrapolation method.

In this subsection, we explore the alternate (exponential) tracer extrapolation method:

$$\begin{split}T(t)=&T(t_{\rm end})\exp(-\frac{t-t_{\rm end}}{3\tau}),\\ \rho(t)=&\rho(t_{\rm end})\exp(-\frac{t-t_{\rm end}}{\tau}),\\ \end{split}$$ (2)

where $t_{\rm end}$ is the post-bounce duration of the 3D simulations and $\tau$ is an expansion timescale fitted using the last 10 ms of the trajectories. The left panel of Figure 7 compares the power-law and exponential methods using an example tracer. The exponential method leads to a much faster decrease of temperature and density than the power-law method. Therefore, nucleosynthesis will terminate significantly earlier if the exponential method is used. We apply the exponential method to the 25 $M_{\odot}$ model, whose power-law version manifests one of the greatest ⁴⁴Ti and ⁵⁷Ni yield growths during the extrapolated phase. The results are shown in the right panel of Figure 7. We find that the choice of extrapolation methods can change the final yield by about 20%. This is not a negligible influence, but a $\sim$20% uncertainty is tolerable in the context of current CCSN nucleosynthesis studies and does not change our general conclusions. The impact of the choice of extrapolation method is weaker for other models (e.g. the 9.x models, 17, 18, and 60 models), because these 3D simulations themselves have already captured much of the late-time yield growth.

A more realistic tracer evolution at late times might be a combination of the power-law and exponential relations. When the tracer is ejected from the surface of the PNS in the neutrino-driven wind, it achieves a radial velocity of a few tens of thousands of kilometers per second. This high velocity leads to a fast temperature and density decrease, with a slope similar to that of the exponential method. When the tracer reaches the wind termination shock where the neutrino-driven wind catches up with the relatively slowly-moving shocked ejecta, its temperature and density evolution switch to the second mode after a short period of shock heating. This second phase has a slope more similar to that obtained using the power-law method. Therefore, an extrapolation method similar to the one employed in Harris et al. (2017) and Sieverding et al. (2023a) might be a better choice, which concatenates the exponential and power-law methods. However, as the position of the wind termination shock is changing, the temperature at which the two methods should be linked is also time-varying. In addition, the wind region can be significantly deformed from a sphere due to interactions with infalling matter (Wang & Burrows, 2023), so there can be a wide range of wind termination shock radii and temperatures, even in a single snapshot of a simulation. Therefore, the only way to unambiguously calculate the late-time evolution of CCSN explosions is to carry the 3D radiation-hydro simulations to 15$-$20 seconds post-bounce.

Figure 7: Comparison of the two thermal history extrapolation methods and their influence on the nucleosynthetic results. Blue lines show the evolution before extrapolation, green lines show the power-law extrapolation, and red lines show the exponential extrapolation. Left: Temperature and density evolution of an example tracer selected from the 25 $M_{\odot}$ model. The temperature of the tracer is about 3 GK at the beginning of the extrapolation. Right: Total ⁴⁴Ti and ⁵⁷Ni produced in the 25 $M_{\odot}$ model as a function of time. It’s clear that the power-law extrapolation allows the tracer to stay at relatively higher temperature for a longer period of time, and thus enhances the yields. Different extrapolation methods lead to about 20% deviation in the final isotopic yields of the 25 $M_{\odot}$ model, which is the model most sensitive to extrapolation choices.

3.3 Stable Isotopes of Nickel and Stable Nickel to Iron Ratios

Not only does the expansion timescale influence the electron fraction ($Y_{e}$) during the early explosion phase, interactions between neutrino-driven winds and fallback accretion can change the outflow velocities and influence the wind $Y_{e}$s. In some models (11, 17, 40 and 60 $M_{\odot}$), we see neutron-rich episodes during the neutrino-driven wind phase. The mass ejected during these periods is very small compared to the mass ejecta during the early explosion phases and the ejecta entropy is higher. These neutron-rich phases can easily produce heavy elements up to ⁹⁰Zr (Wang & Burrows, 2023, 2024a). One observable in supernova remnants, the stable Ni/Fe ratio, is dependent upon the strength of such neutron-rich phases. There are roughly three classes of models in our set. The first class contains the most neutron-rich models in which the early explosion phase has $Y_{e}<0.5$ (e.g., model 9a, and the models in Wang & Burrows (2024b)). These models have Ni/Fe ratios more than 7 times as high as the solar value, and are close to the recent observation of the Crab nebula (Temim et al., 2024). The second class contains the models with neutron-rich episodes in their winds (11, 17, 40, and 60 $M_{\odot}$). The Ni/Fe ratios in these models are about 1.5 to 5 times as high as the solar value, depending upon the amount of mass ejected during the neutron-rich period. The last class represents all models without any neutron-rich period. These models all have sub-solar Ni/Fe ratios (as low as 0.6 times the solar value). Although the neutron-rich wind phases play an important role in the production of many interesting isotopes, their appearance depends upon many aspects of the explosion dynamics and is very difficult to model or discern in advance of simulation. We list the asymptotic ratios of stable nickel to stable iron that we find for this suite of 18 long-term 3D CCSN models in Table 1. The shape and direction of the wind regions and accretion funnels, the PNS properties, and the outer envelope structures may all have an impact on the wind electron fraction. Therefore, we can conclude only that the neutrino-driven winds are mostly proton-rich, but with episodic neutron-rich episodes.

4 Conclusion

This paper presents the late-time evolution of the yields of ⁴⁴Ti, ⁵⁶Ni, ⁵⁷Ni, and stable iron and nickel isotopes for 18 3D CCSN simulations. The ejecta isotopic abundances at about twenty seconds post-bounce are calculated, the last phases by an extrapolation method of both tracer thermal histories and neutrino-driven wind mass outflow rates from our already long-term 3D simulations. For the first time, a detailed systematic study reveals the theoretical asymptotic yields of these interesting nucleosynthetic products, informed by state-of-the-art long-term 3D CCSN simulations.

In our models, the majority of ⁵⁶Ni is produced in the first 4$-$5 seconds after bounce, revealing that the explosive nucleosynthetic component and the early-phase neutrino-driven wind dominate ⁵⁶Ni production. A clear relation between final ⁵⁶Ni yield and explosion energy has previously been shown in Wang & Burrows (2024a) and Burrows et al. (2024).

We find that the production timescales of ⁴⁴Ti are significantly longer than those of ⁵⁶Ni. Some models take more than 15 seconds to reach their asymptotic ⁴⁴Ti abundances. We find that ⁴⁴Ti is mostly produced in neutrino-driven winds via the “(p,$\gamma$)-leakage” channel (Magkotsios, 2011). In this special proton-rich variant of freeze-out, the (p,$\gamma$) reactions transfer the nucleosynthetic flow from symmetric to proton-rich isotopes. As the temperature decreases, the freeze-out process starts and weak interactions decrease $Y_{e}$ and cause the flow to revert to more stable isotopes. In this process, isotopes such as ⁴⁵V and ⁴⁶Cr are produced at early times, and they are slowly converted back to ⁴⁴Ti as the temperature decreases. This freeze-out process remains active even when the temperature is lower than 1 GK. Therefore, there are two main reasons for the long production timescale of ⁴⁴Ti. First, the neutrino-driven wind phase in which ⁴⁴Ti is produced lasts for many seconds post-bounce. The simultaneous explosion and accretion phenomenon in 3D CCSN simulations allows the outer stellar matter to fall very close to the PNS through accretion funnels, which significantly increases the amount of mass interacting with neutrinos and leads to extended and stronger neutrino-driven winds (Wang & Burrows, 2023). This enhancement of wind strength lasts for many seconds until accretion is completely terminated and a spherical wind develops (Stockinger et al., 2020; Wang & Burrows, 2024b). Second, the freeze-out process requires a very low temperature (below 1 GK) to reach the asymptotic abundances, which also takes many seconds to achieve.

Our 3D models produce ⁴⁴Ti ranging from $10^{-6}M_{\odot}$ to $2\times 10^{-4}M_{\odot}$. The high ⁴⁴Ti/⁵⁶Ni ratio observed in Cas A can easily be achieved by many of our models. Therefore, the “⁴⁴Ti-problem” (The et al., 2006) is no longer an issue in self-consistent long-term, 3D, initially non-rotating CCSN models.

We find that there are various pathways to ⁵⁷Ni production. If the explosion has ever experienced a neutron-rich phase, the ⁵⁷Ni/⁵⁶Ni ratio quickly rises to a high value during the neutron-rich phase and stays there or decreases slowly due to ⁵⁶Ni production in the neutrino-driven wind. For models with no neutron-rich ejecta, the production timescale of ⁵⁷Ni is even longer than that of ⁴⁴Ti. This is because ⁵⁷Ni is also produced in significant amounts in the $\alpha$-rich freeze-out process in proton-rich neutrino-driven winds. The asymptotic ⁵⁷Ni/⁵⁶Ni ratios for the two pathways are not as different as they are initially, unless the neutron-rich phase is early and strong.

The observed ⁴⁴Ti/⁵⁶Ni and ⁵⁷Ni/⁵⁶Ni ratios in SN1987A are both $\sim$2 sigma above the values seen in our 3D simulations. One possible explanation is that the actual explosion might have an early and strong neutron-rich phase, so that a fraction of the explosive nucleosynthetic products become neutron-rich. This will enhance ⁵⁷Ni production, but the most significant influence is that the production of ⁵⁶Ni at the early phase will be suppressed. ⁴⁴Ti is not influenced by this early neutron-rich phase, since it is produced mostly in the wind phase. Therefore, a more neutron-rich explosion, which might be caused by a more prompt explosion, will naturally lead to higher ⁴⁴Ti/⁵⁶Ni and ⁵⁷Ni/⁵⁶Ni ratios.

To test the robustness of the above findings, we changed our tracer extrapolation method from a power-law method to an exponential method. We find that the choice of extrapolation method changes the final yields by no more than 20%, even in our worst case. This is not a negligible influence, but a $\sim$20% uncertainty is tolerable for current CCSN nucleosynthesis studies and it will not change our general and substantive conclusions. However, the only way robustly to determine the late-time nucleosynthetic yields of core-collapse supernovae is to carry out detailed 3D supernova simulations to 15$-$20 seconds post-bounce. This project, though daunting, is now almost within reach of modern codes and machines, but will require a concerted effort to realize.

We also briefly discuss the stable Ni/Fe ratios. Our models have stable Ni/Fe ratios for all the ejecta of the explosion ranging from sub-solar values (as low as 0.6 times the solar value) to seven times the solar values. The latter is similar to the recent observation value of the Crab nebula (Temim et al., 2024). However, Ni/Fe ratios, and in general the production of neutron-rich isotopes, depend strongly upon the strength of neutron-rich phases in the explosions. Since such phases are sensitive to many aspects of explosion dynamics, their occurrence seems to be stochastic. Therefore, the uncertainties of neutron-rich isotope abundances can be very large.

It is important to list some limitations of this work. First, our results rely on extrapolations of tracer thermal histories and neutrino-driven wind outflow rates. The accuracy of these extrapolations has not been tested by 3D simulations carried out to the same timescale. Second, although this is a systematic study covering a number of different progenitor models, all progenitors are calculated using the code Kepler (Sukhbold et al., 2016, 2018) and it is unclear if these progenitors cover the same range of structures and properties of real stars. In addition, the nuclear reaction rates are all taken from the JINA Reaclib (Cyburt et al., 2010), for which there are still large uncertainties in some rates (Bliss et al., 2020). Although we have mentioned briefly the potential roles of initial perturbations and rotation, a detailed study of their effects is still lacking. There are also uncertainties in the nuclear equation of state and the neutrino opacities/emissivities, etc.

Data Availability

This paper uses previously published data in https://doi.org/10.5281/zenodo.10498614. This online repository includes all the tracer trajectories (masses, positions, temperatures, electron fractions, and local neutrino spectra), together with the ejecta abundances of each model at the end of the 3D simulations. The progenitor models from Sukhbold et al. (2016) and Sukhbold et al. (2018) are available as well. The 17, 18, 18.5, 20, 25, and 60 $M_{\odot}$ models have been evolved to later times in 3D since the publication of the above repository. However, we can not replace the files in the repository, as this would significantly exceed to size limit of the repository. Therefore, we encourage readers to start with the online repository, and if necessary, to make direct requests to the authors for the rest of data.