Nucleosynthesis for SN 1987A from single-star and binary-merger progenitors

SN 1987A is the best-studied supernova. It is a core-collapse supernova (CCSN) that occurred in the Large Magellanic Cloud and was first detected on Earth on February 23.316 UT, 1987. Even before the first light had arrived, an initial neutrino burst was detected by Kamikande II [1, 2], the IMB-detector [3, 4], and possibly also by the Baksan Neutrino Observatory [5, 6]. The early optical light curve brightened very rapidly (by a factor of 100 over ∼3 h) indicating a more compact atmosphere than is typical for red supergiants (RSGs). Early spectra revealed hydrogen features at high velocities (up to 0.1c), classifying it as a Type II supernova.

Many groups have used single-star evolutionary models to compute the ejecta mass and velocities as a function of the explosion energy and to fit the observed light curve of SN 1987A [7–15]. From such studies, estimates for the explosion energy, progenitor mass, and ejecta masses of ⁵⁶Ni are obtained. For example, Blinnikov et al [11] determined an explosion energy of ${E}_{\mathrm{expl}}=(1.1\pm 0.3)\times {10}^{51}$ erg assuming ∼14.7 M_⊙ of ejecta and 10.3 M_⊙ of hydrogen-rich envelope. Seitenzahl et al [15] determined the yields of ⁵⁶Ni and ⁵⁷Ni from a least-squares fit of the decay chains to the observed bolometric light curve. The amount of ⁴⁴Ti ejecta can be estimated from the energy deposition by positrons originating from the decay of ⁴⁴Sc (the daughter nucleus of ⁴⁴Ti) [15]. However, observation of the 67.9 and 78.4 keV x-ray lines from the decay of ⁴⁴Ti with INTEGRAL [16] or NUSTAR [17] is a more direct approach. Boggs et al inferred $(1.5\pm 0.3)\times {10}^{-4}\,{M}_{\odot }$ of ⁴⁴Ti from the NuSTAR data.

From pre-explosion imaging, the progenitor star of SN 1987A has been identified as the compact B3 supergiant Sanduleak −69° 202 (Sk -69 202) [18–21]. The deduced surface temperature is ${\mathrm{log}}_{10}{T}_{\mathrm{eff}}\approx 4.2$ [7], which corresponds to a blue supergiant (BSG) at the time of explosion. However, the evolutionary pathway to the BSG stage is still not well understood. Stellar evolution theory predicts that massive stars end their lives as RSGs, unless exotic scenarios are invoked or the parameters are finely tuned (see [9, 22, 23] for reviews). An alternative approach is to consider stellar evolution in binary systems. The polarization of light from SN 1987A indicates a flattened envelope structure, which is not consistent with a single-star envelope [24]. The triple-ring structure of the nebula of SN 1987A is also difficult to explain with single-star evolution. Moreover, the inner ring is enhanced in CNO-processed material, which is also not expected from single-star evolution. A merger scenario in a binary system could explain the BSG nature of Sk -69 202 prior to explosion [25, 26]. The best current understanding of such a scenario involves the complete merger of the stellar components of a wide binary system, in which the primary star fills its Roche lobe only after core helium burning is completed (see [27]). Recently, two groups created stellar evolution models for the progenitor of SN 1987A based on this merger scenario [28, 29], which fulfill the observational constraints from Sk -69 202 and which are now available for use in explosion and nucleosynthesis simulations of SN 1987A.

The modeling of CCSNe traditionally follows one of two approaches: Sophisticated, computationally very expensive multi-dimensional simulations are used to study the details of the explosion mechanism. Artificially induced explosions based on a piston or thermal bomb approach are employed to predict nucleosynthesis yields. In recent years, several groups have developed effective explosion models for CCNSe as complimentary approach [30–33]. These effective models mimic in spherical symmetry the net effects of realistic, three-dimensional simulations to trigger explosions in otherwise non-exploding simulations that remain computationally efficient to allow for large numbers of simulations. In particular, they follow the evolution from pre-explosion through collapse, bounce, and onset of explosion, and they ensure an accurate treatment of the electron fraction by including a spectral neutrino transport scheme and following the detailed evolution of the proto-neutron star. This is of particular importance for the innermost ejecta where the electron fraction can be significantly altered from its pre-explosion value, and hence also the resulting nucleosynthesis [34–36]. Such an approach is suitable for questions that require many simulations, for example investigating the connection between the progenitor and the explosion outcome, successful or failed explosion, or predicting accurate nucleosynthesis yields for the entire mass range of CCSNe. Being effective models, as opposed to fully self-consistent, there are parameters that require calibration using observational data. Such effective and semi-analytic models have been used to study explosion properties and outcomes [30–33, 37, 38]. In two cases, also detailed nucleosynthesis yields were computed [36, 39]. In this work, we combine our effective CCSN framework, the PUSH method [31, 33], with the recent simulations for binary progenitors of SN 1987A [28] to compute nucleosynthesis yields for SN 1987A and to compare to our earlier results obtained for core-collapse progenitors from single-star evolution.

The paper is organized as follows: in section 2, we describe the pre-explosion models and our computational setup. We also summarize the key aspects of the PUSH method used to trigger explosions. The calibration procedure is explained in section 3. We present the detailed nucleosynthesis yields in section 4 and compare the results to those from single-star pre-explosion models. We conclude with a discussion of our results in section 5.

2.1. Initial models

For this study, we use two sets of stellar evolution calculations, MH17 [28] and WHW02 [40], as starting points for our CCSN explosion and nucleosynthesis simulations. We refer to these starting points as 'pre-explosion models'. A list of all the pre-explosion models discussed in this work is given in table 1. Each pre-explosion model has a unique label, for which we follow the same naming conventions as in the original publications. For the MH17 series, the letter indicates the amount of mass of the He core that is dredged up: 'a' corresponds to 10% of He-shell mass, 'b' to 50% of the He-shell mass. For the WHW02 series, the letter 's' indicates solar metallicity. For both series, the number(s) encode(s) the ZAMS mass(es).

Table 1.  Pre-explosion models considered in this study.

Series	M1	M2	Final mass	Compactness ξ2.0	Label	Reference
	( M⊙)	( M⊙)	( M⊙)	(−)
MH17	16.0	4.0	19.0	0.298 43	a16-4	[28]
MH17	17.0	7.0	22.8	0.289 12	a17-7	[28]
MH17	17.0	8.0	23.8	0.439 02	a17-8	[28]
MH17	15.0	7.0	21.1	0.221 54	b15-7	[28]
MH17	15.0	8.0	22.1	0.124 99	b15-8	[28]
MH17	16.0	7.0	22.0	0.172 49	b16-7	[28]
WHW02	18.8	—	15.0	0.244 05	s18.8	[40]

One set, MH17, represents the outcome of the binary merger of a primary RSG star and a secondary main-sequence (MS) star. As the primary star expands to a RSG at the end of core helium burning, it fills and overflows its Roche lobe. This results in unstable mass transfer, so-called 'mass transfer of Type C' [41], and eventually the secondary MS star is completely dissolved in the common envelope of the primary RSG. The material of the secondary MS star is mixed with the convective envelope of the primary RSG star. Some H-rich material of the secondary star can even penetrate the He core of the primary RSG and dredge up nitrogen-rich and helium-rich material. Mass transfer ends when the secondary star is completely dissolved in the envelope of the primary RSG. This results in stellar structures with some enrichment of the envelope with material from the He core and usually a larger ratio of envelope-mass to core-mass than is seen in single stars. The entire MH17 set consists of 84 pre-explosion models, of which 59 represent hot BSGs at the end of their evolution. Here, we include only the six pre-explosion models from MH17 that are hot BSGs and also fulfill the observational constraints of SK -69 202, in particular the temperature, luminosity and surface abundances (table 4 in [28]).

The other set, WHW02, are spherically symmetric, non-rotating, solar metallicity pre-explosion models, computed with the stellar evolution code Kepler. This set consists of 101 individual pre-explosion models with ZAMS masses from 10.8 M_⊙ to 75 M_⊙, which represent the evolution of single stars with evolutionary end points of RSGs. Pre-explosion models from WHW02 are readily available to the community, and hence, they are often used in CCSN simulations, e.g. [42–50]. We used a subset of WHW02 (their models with ZAMS masses of 18–21 M_⊙) to represent the progenitor of SN 1987A in previous work [33]. Here, we rely on the best-fit model (s18.8) which we found in our earlier work [33].

Figure 1 shows the density profiles of all the pre-explosion models from table 1. The largest differences are in the envelope, which is not surprising given that MH17 are BSGs whereas WHW02 are RSGs. In the innermost 3 M_⊙, however, all of them are quite similar. It is this innermost region that affects the explosion properties and the nickel yields the most. A useful quantity to characterize the structure of the pre-explosion models in a single number is the compactness [51], which describes the ratio of a given mass, M, and the radius, R(M), which encloses that mass:

$$\begin{eqnarray}&&{\xi }_{M}\equiv \displaystyle \frac{M/{M}_{\odot }}{R(M)/1000\,\mathrm{km}},\end{eqnarray} \tag{ 1 }$$

where we use $M=2.0\,{M}_{\odot }$ as before. The pre-explosion models in table 1 have a similar compactness values with the exception of a17-8 (higher compactness) and b15-8 (lower compactness).

Zoom In Zoom Out Reset image size

2.2. Explosion simulations

For the explosion simulations, we use Agile, a general-relativistic, adaptive-mesh hydrodynamics code, which is coupled to spectral neutrino transport (IDSA for electron-flavor neutrinos [52] and ASL for heavy-flavor neutrinos [53]). For the nuclear equation of state we use HS(DD2) for matter in NSE conditions [54] and an ideal gas approach coupled to an approximate alpha-network for non-NSE conditions [31]. We perform our simulations in spherical symmetry using the PUSH method [31, 33, 36]. The PUSH method triggers explosions in spherical symmetry by artificially enhancing the neutrino heating. This is achieved by depositing a fraction of the heavy-lepton flavor neutrinos in the heating region behind the shock via a local heating term

$$\begin{eqnarray}&&{Q}_{\mathrm{push}}^{+}(t,r)=4\,{ \mathcal G }(t){\int }_{0}^{\infty }{q}_{\mathrm{push}}^{+}(r,E)\,{\rm{d}}{E},\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray}&&{q}_{\mathrm{push}}^{+}(r,E)\equiv {\sigma }_{0}\,\displaystyle \frac{1}{4\,{m}_{b}}{\left(\displaystyle \frac{E}{{m}_{e}{c}^{2}}\right)}^{2}\displaystyle \frac{1}{4\pi {r}^{2}}\left(\displaystyle \frac{{{\rm{d}}{L}}_{{\nu }_{x}}}{{\rm{d}}{E}}\right){ \mathcal F }(r,E).\end{eqnarray} \tag{ 3 }$$

This heating depends on the spectral energy flux $({{\rm{d}}{L}}_{{\nu }_{x}}/{\rm{d}}{E})/(4\pi {r}^{2})$ for any single flavor of neutrinos, ν_x , with energy E. The other terms are the typical neutrino cross-section (${\sigma }_{0}\approx 1.759\times {10}^{44}$ cm²) and the average baryon mass (${m}_{b}\approx 1.674\times {10}^{-24}\,{\rm{g}}$). The spatial dependence of PUSH is in the function ${ \mathcal F }(r,E)$, which ensures that the additional heating only takes place where electron neutrinos are also heating. The temporal evolution of PUSH is in the function ${ \mathcal G }(t)$. This function ${ \mathcal G }(t)$ is zero until time t_on = 0.8 s post-bounce, when it increases linearly from zero to its maximum value ${k}_{\mathrm{push}}$  over the time interval t_rise. Then, ${ \mathcal G }(t)$ remains constant at ${k}_{\mathrm{push}}$ until, at 1.0 s post-bounce, it decreases linearly from ${k}_{\mathrm{push}}$ back to zero over the same time interval t_rise. The two free parameters in the method are ${k}_{\mathrm{push}}$  and ${t}_{\mathrm{rise}}$. This is exactly the same setup as used in [33, 36]. We refer the reader to these papers for a more detailed discussion.

Our simulations include all the matter of the pre-explosion model from the center to the He-layer, only excluding the low-density H-envelope. We follow the collapse, the formation of the proto-neutron star (PNS), and the subsequent explosion. We run the simulations for a total simulation time of 5 s, corresponding to 4.6 s–4.8 s post bounce. In this setup, the mass, i.e. the bifurcation between matter ejected and matter settling back onto the PNS, emerges naturally from the simulation. It is identified following the suggestion by Bruenn in [55] as the mass outside the homologous core, which maximizes the explosion energy (see also equation (17) in [31]).

In the innermost ejecta, the details of the nucleosynthesis depend strongly on the electron fraction (Y_e). Here, weak processes during collapse, bounce, and explosion significantly change the Y_e from its pre-explosion value. Of particular importance are the charged-current neutrino reactions ${\rm{n}}+{\nu }_{e}\leftrightarrow {\rm{p}}+{{\rm{e}}}^{-}$ and ${\rm{p}}+{\bar{\nu }}_{e}\leftrightarrow {\rm{n}}+{{\rm{e}}}^{+}$ [56–58]. In multi-dimensional simulations, small changes in the neutrino-matter coupling, for example due to many-body corrections in the neutrino-nucleon cross sections, can amplify near critical conditions [59]. However, in spherically symmetric simulations this sensitivity has not been seen. Similar to multi-dimensional simulations, we do not include neutrino flavor transformations in our simulations.

2.3. Nucleosynthesis calculations

The detailed nucleosynthesis yields are calculated using a postprocessing approach with the full nuclear reaction network CFNET. For the reaction rates, we use the REACLIB library [60] which includes experimentally known rates whenever available and statistical model predictions for n-, p-, and alpha-capture reactions otherwise [61]. For weak interactions we include electron and positron capture reactions [62], β⁻and β⁺ decays from the nuclear database NuDat2 (where available) and theoretical predictions from Möller et al [63] otherwise. Neutrino and anti-neutrino captures on free nucleons are also included [35].

We post-process the ejecta in mass elements ('tracers') of ${10}^{-3}\,{M}_{\odot }$ each. The hydrodynamic history of each tracer is obtained from the hydrodynamics simulation. If the temperature of the tracer at the end of the hydrodynamic simulation is still high enough for nucleosynthesis to occur, we extrapolate the hydrodynamic history using a free expansion for the density and adiabatic expansion for the temperature [64]:

$$\begin{eqnarray}&&r(t)={r}_{\mathrm{final}}+{{tv}}_{\mathrm{final}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\rho (t)={\rho }_{\mathrm{final}}{\left(\displaystyle \frac{t}{{t}_{\mathrm{final}}}\right)}^{-3},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&T(t)=T[{s}_{\mathrm{final}},\rho (t),{Y}_{e}(t)].\end{eqnarray} \tag{ 6 }$$

In these equations, r is the radial position, v is the radial velocity, ρ is the density, T is the temperature, s is the entropy per baryon, and Y_e is the electron fraction of the mass zone. The subscript 'final' refers to the end of the hydrodynamic simulation. The temperature is calculated using the Timmes & Swesty EOS [65].

In this section, we discuss how the free parameters in the PUSH method are calibrated using observational data from SN 1987A. The general idea is to find a combination of a pre-explosion model (with ZAMS mass matching the progenitor mass of SN 1987A) and values for the free parameters ${k}_{\mathrm{push}}$ and ${t}_{\mathrm{rise}}$ such that all observational quantities of SN 1987A are simultaneously reproduced.

In the present work, we adopt the same observationally derived quantities for SN 1987A as in [33]: for the explosion energy we use ${E}_{\mathrm{expl}}=(1.1\pm 0.3)\times {10}^{51}$ erg as reported in [11]. A detailed summary of explosion energy estimates can be found for example in [66]. The adopted estimate was obtained based on the assumption of ∼14.7M_⊙ of ejecta and ∼10.3 M_⊙ of hydrogen-rich envelope. For the element abundances of ⁵⁶Ni and ⁵⁷Ni, we adopt the values obtained from a least squares fit to the bolometric light curve [15]. The yield of ⁵⁸Ni was taken from [67]. Different techniques have been used to estimate the yield of ⁴⁴Ti [15–17, 67–70]. Here, we adopt a yield of $(1.5\pm 0.3)\times {10}^{-4}\,{M}_{\odot }$ obtained from NUSTAR observations [17]. The mass of the progenitor of SN 1987A has been determined to be in the range of 18M_⊙–21 M_⊙ [7, 10].

Below we will discuss the calibration procedure using the example of the MH17 set. The same calibration procedure was already preformed in earlier work [31, 33] for the WHW02 set and resulted in finding a best-fit model for SN 1987A within the WHW02 set which we use in this work for the comparison presented in section 4.

For the calibration procedure, we performed hydrodynamic simulations for each of the six pre-explosion models of the MH17 set, using the setup described in section 2.2 for a range of values for the free parameters ${k}_{\mathrm{push}}$  and ${t}_{\mathrm{rise}}$. For each simulation, we computed the explosion energy, assuming that the rest-mass-subtracted total energy of the ejecta is eventually converted into kinetic energy, and we computed the detailed nucleosynthesis yields using the nuclear reaction network. A subset of all the simulations performed is shown in figure 2. Simulations which resulted in explosion energies much lower, i.e. to the left of the error box, or much higher, i.e. to the right of the error box, than the adopted range for SN 1987A (shown as error box) are omitted from the figure for clarity.

Zoom In Zoom Out Reset image size

One can observe several general trends from these results: when looking at the ⁵⁶Ni yields (top left panel of figure 2), we find that the six pre-explosion models fall into three groups: only for three models (a16-4 in red, a17-7 in green, and b15-7 in blue) can we find combinations of ${k}_{\mathrm{push}}$  and ${t}_{\mathrm{rise}}$  that result in an explosion energy and a ⁵⁶Ni yield within the observed ranges. When requiring explosion energies consistent with observations, one model (a17-8 in purple) always results in ⁵⁶Ni yields that are too high and two models (b15-8 in teal and b16-7 in brown) give ⁵⁶Ni yields that are too low. For each of the models (see all points in a single color), one notices a trend that lower explosion energies correlate with lower ⁵⁶Ni yields. We previously found the same correlation within the WHW02 set [31, 33]. Due to this correlation, it is sufficient to investigate simulations from combinations of ${k}_{\mathrm{push}}$  and ${t}_{\mathrm{rise}}$  that result in explosion energies within the required range of (1.1 ± 0.3 ) × 10⁵¹ erg. As seen in the other panels of figure 2, models that do not match observations in ⁵⁶Ni yields also do not have yields matching observations for the other three isotopes. For the two asymmetric Ni isotopes (⁵⁷Ni in the top right panel, ⁵⁸Ni in the bottom left panel), the yields are lower than the observed values. The only exception is model b15-7 (blue) which yields ⁵⁷Ni within the observed range and ⁵⁸Ni above the observed value. However, note that for ⁵⁸Ni the observed value is quite uncertain and no error bar has been reported [67]. In all our simulations, the ⁴⁴Ti yield (bottom right panel of figure 2) is much lower than the observed value. As such, the observed ⁴⁴Ti yield does not provide any constraint on the calibration procedure, however we include it for completeness as observational yields are available only for the four isotopes discussed (^56,57,58Ni and ⁴⁴Ti). From this analysis, we conclude that only b15-7 can be calibrated within our framework to represent SN 1987A. We adopt b15-7 together with k_push = 4.1 and t_rise = 500 ms as our model for SN 1987A from the MH17 set.

In previous work [33, 36] we performed the same calibration procedure as described above for the WHW02 pre-explosion models with ZAMS masses of 18–21 M_⊙. The trends and general behavior were the same as discussed above for the MH17 set (we refer the interested reader to [31, 33] for details). Our analysis resulted in adopting model s18.8 and values of k_push = 4.3 and t_rise = 400 ms as the calibration model for SN 1987A, which we use in this work without repeating the calibration procedure.

We summarize the parameters of our SN 1987A models in table 2 and show the resulting yields for ^56,57,58Ni and ⁴⁴Ti together with their observed values in figure 3. In addition to the two best-fit models, figure 3 also shows simulations for b15-7 with parameter values of ${k}_{\mathrm{push}}$  and ${t}_{\mathrm{rise}}$ which are slightly different from their best fit values (±10 ms for ${t}_{\mathrm{rise}}$  and ±0.1 for ${k}_{\mathrm{push}}$). From these results, one can see that the ⁵⁶Ni yield represents the strongest constraint on the calibration for a given range of explosion energies. Moreover, one can see that the nucleosynthesis yields can distinguish between different models of similar explosion energy.

Zoom In Zoom Out Reset image size

Table 2.  Best fit parameters for all models able to reproduce SN 1987A (see text for details). The first two lines summarize the observational quantities^a of SN 1987A.

Name	kPUSH	trise	Eexpl	$m{(}^{56}\mathrm{Ni})$	$m{(}^{57}\mathrm{Ni})$	$m{(}^{58}\mathrm{Ni})$	$m{(}^{44}\mathrm{Ti})$
	(–)	(ms)	(1051erg)	( M⊙)	( M⊙)	( M⊙)	( M⊙)
SN 1987A	—	—	1.1	0.071	0.004 1	0.006	1.5 × 10−4
	—	—	±0.3	±0.003	±0.001 8		±0.3 × 10−4
b15-7	4.1	500	1.1	0.072	0.003 0	0.007 4	3.60 × 10−5
s18.8	4.3	400	1.2	0.069	0.002 7	0.006 6	3.05 × 10−5

Note.

^aE_expl: [11]; ^56,57Ni: [15]; ⁵⁸Ni: [67]; ⁴⁴Ti: [17].

In the previous section we have discussed how we have identified the two pre-explosion models that can represent SN 1987A and the parameter values for ${k}_{\mathrm{push}}$ and ${t}_{\mathrm{rise}}$ (see table 2) for each of them. Now we turn our attention to the detailed nucleosynthesis predictions from these two simulations. Figure 4 shows the integrated final abundances for both models. It is noteworthy that only the iron group yields are strongly affected by the explosion. The intermediate and lighter elements are more strongly dependent on the composition of the pre-explosion model.

Zoom In Zoom Out Reset image size

Iron group nuclei are synthesized in layers that undergo either complete or incomplete silicon burning. The nucleosynthesis in these layers is sensitive to the location of the mass cut which determines how much of the silicon–oxygen shell is ejected. The nucleosynthesis is also sensitive to the weak interactions such as electron capture reactions and neutrino-induced reactions during the collapse and the explosion phase, as these reaction alter the electron fraction from its pre-explosion value. In the innermost layers this change can be quite significant, from slightly neutron-rich conditions in the pre-explosion state to proton-rich conditions in the post-explosion phase. The nickel isotopes ^56,57,58Ni (as well as ⁴⁴Ti) are produced in the innermost ∼1.2 M_⊙ of ejecta in our simulations. The total yields of ⁵⁶Ni are very similar between b15-7 and s18.8 due to our calibration method (see also figure 3) which requires (0.071 ± 0.003) M_⊙ of ⁵⁶Ni. Moreover, b15-7 and s18.8 have similar compactness values ξ_2.0 at bounce, so it is not surprising that for similar explosion energies the ⁵⁶Ni yields are also similar. A correlation between compactness and Ni yield has been discussed in [31, 33]. The mass cut is in the silicon–oxygen layer at Y_e of 0.497 9 in both cases (figure 5) and ⁵⁶Ni is synthesized through explosive burning of Si. In the case of b15-7, ⁵⁶Ni is also co-produced—at the few percent level—with Si further outwards in the star. The resulting extended Ni distribution may impact the light curve through the so-called 'Ni-bubble effect' [71, 72]. The yields of isotopes with $56\leqslant A\leqslant 60$, which are co-produced in the same layers, are also slightly higher for b15-7 compared to s18.8. In order for the neutron-rich nickel isotopes ^57,58Ni to be synthesized in sufficient amounts, slightly neutron-rich conditions (Y_e < 0.5) are required such that alpha-rich freeze-out results in nuclei one to two nucleons away from the N = Z line, as is the case for b15-7 and s18.8. The yield of ⁴⁴Ti is very sensitive to the location of the mass cut [73], which in the PUSH framework emerges from the simulation and is hence not an additional free parameter (unlike in traditional piston or thermal bomb models). In both models, the ⁴⁴Ti yield is about a factor of 5 lower than the observed value. In previous work [31], we estimated the increase in ejected ⁴⁴Ti in our models by assuming homogeneous mixing in the layers between the mass cut and the outer boundary of the Si-rich shell and assuming some fallback by hand. See [74] for the mixing-fallback model. We found that an increase by a factor of 2–3 is possible, which is considerable but not sufficient to match observations.

Zoom In Zoom Out Reset image size

For nuclei lighter than iron, the largest differences are in the mass range of 20 ≲ A ≲ 40, where s18.8 resulted in higher yields for almost all isotopes when compared to b15-7. These nuclei are synthesized in layers that reach peak temperatures of less than 5 GK, but still high enough for explosive processing of the O-rich pre-explosion material (figure 5). For b15-7, most of the ejected Si originates from the oxygen-rich pre-explosion layer and is freshly synthesized during the explosion. In this case, the amount of Si ejecta is limited by how much of the O-rich pre-explosion material reaches high enough temperatures during the explosion for explosive oxygen burning. In addition, the abundance profiles of Si (and S) are extended further outwards in the star for b15-7 at the ∼6% (∼0.5%) level. In the case of s18.8, however, the ejected Si has been made prior to the explosion and is ejected unprocessed. As a result, s18.8—which has a thick pre-explosive Si-layer—has higher Si abundances in the ejecta.

The total yields of ¹⁶O and ¹²C yield are higher for b15-7 than for s18.8. Both isotopes originate from the O-rich layer. Whereas this O-rich layer is narrower for b15-7 than for s18.8, the mass fractions of oxygen (∼83% for b15-7 compared to ∼75% for s18.8) and also of carbon are higher for b15-7, resulting in overall higher final yields of ¹⁶O and ¹²C.

Figure 4 also reveals that some amount of nuclei heavier than iron are synthesized. This happens in the innermost ejecta. We find proton-rich conditions (Y_e ∼ 0.52–0.53) in 0.6–0.8 M_⊙ of inner ejecta of b15-7 and s18.8. In these proton-rich layers, the composition is dominated by ⁵⁶Ni, ⁴He, and protons (as expected for proton-rich, alpha-rich freeze-out) together with ⁶⁰Ni and ⁶⁴Zn, which are produced chiefly as ⁶⁰Zn and ⁶⁴Ge. Under slightly neutron-rich conditions (Y_e < 0.5) isotopes with A > 80 are produced. We find conditions of Y_e < 0.5 in the very innermost ejected tracers. However, our mass resolution of 10⁻³ M_⊙ per tracer is too coarse to accurately resolve this late-time neutrino-driven ejecta. The abundances for A > 80 seen in figure 4 originate from a few tracers with Y_e ∼ 0.43, outside of which the Y_e steeply rises from ∼0.43 to ∼0.52 over only one or two tracers, corresponding to $(1-2)\times {10}^{-3}\,{M}_{\odot }$. If this tracer happens to be placed at Y_e ∼ 0.45, this may result in considerable contributions to ⁶⁶Zn.

Finally, we compare the iron group yields from the two SN 1987A models (red circles for s18.8; blue stars for b15-7) to the observationally derived abundances for the metal-poor star HD 84937 (figure 6). These abundances have been derived using improved laboratory data for iron group transitions [75]. The pattern is quite similar for both models (except for zinc) and both models are overall in good agreement with the observational data. Zinc is co-produced with iron group nuclei as ⁶⁴Zn in proton-rich layers. In layers with ${Y}_{e}\approx 0.45$ zinc can also be produced as ⁶⁶Zn. Figure 6 also includes simulations for b15-7 with slightly different parameter values than the best fit (grey symbols; see also figure 3).

Zoom In Zoom Out Reset image size

In this work, we have used the PUSH framework to model the detailed nucleosynthesis of SN 1987A. For this, we used a set of six pre-explosion models (MH17) representing BSGs at the onset of collapse and matching the observed properties SK-69 202 (the progenitors of SN 1987A). We calibrated the free parameters of the PUSH method by requiring that the simulations simultaneously reproduce the explosion energy and the observed yields of ^56,57,58Ni and ⁴⁴Ti. Of the six models in MH17, this was only possible for b15-7.

For the other five models it was not possible to find any parameter values that reproduced the observed quantities. The model with the highest compactness (a17-8; ξ_2.0 = 0.43902) yielded between 0.109 and 0.112 M_⊙ of ⁵⁶Ni, which is consistent with the expectation of 0.1 M_⊙ from the correlation between ⁵⁶Ni yield and compactness we found in earlier work using the PUSH method, but which is not consistent with the observed value. Within the PUSH framework, it is not possible to get significantly lower ⁵⁶Ni yields without requiring additional late-time fallback by hand for models with such a high compactness. Models b15-8 and b16-7, which have very narrow Si-layers (∼0.01 M_⊙), synthesize most of the ejected ⁵⁶Ni from explosive oxygen burning, resulting in significantly less ⁵⁶Ni than the observationally required (0.071 ± 0.003) M_⊙. Moreover, the mass cut always lies around the interface between the Si-rich and the O-rich layer, which further reduces the amount of ejected Ni. For models a16-4 and a17-7, we found parameters that give the desired explosion energy and ⁵⁶Ni; however, the yields of the asymmetric Ni isotopes ⁵⁷Ni and 58Ni are too low. For these two models the mass cut always resides in the O-rich layer, where the pre-explosion Y_e is 0.4999. This is higher than the Y_e in the Si-rich layer (Y_e ≈ 0.4979) where the mass cut typically is located. Due to this it is not surprising that the yields of the non-symmetric Ni isotopes (⁵⁷Ni and ⁵⁸Ni) are lower, as they require slightly neutron-rich conditions to be synthesized. In particular, ⁵⁸Ni is lower by a factor of 10 compared to the best fit models. In return, isotopes with mass number $A\geqslant 64$ are enhanced, however they only account for a small fraction of the mass in these layers.

In previous work, we performed the same calibration procedure with the WHW02 pre-explosion model set, which is based on single-star evolution at solar metallicity and which represents RSGs at the onset of collapse. There, we found that s1.8.8 with k_push = 4.3 and t_rise = 400 ms is able to simultaneously reproduce the observed explosion energy and ^56,57,58Ni yields of SN 1987A.

We compared the nucleosynthesis yields obtained for b15-7 (k_push = 4.1 and t_rise = 500 ms) with those obtained for s18.8. They both have similar compactness ξ_2.0. As a results, the values of the free parameters ${k}_{\mathrm{push}}$ and ${t}_{\mathrm{rise}}$ are also very similar. Overall, the nucleosynthesis yields, including the amount of ⁵⁶Ni, are very similar. In the iron group, the only noteworthy difference is found for zinc where [Zn/Fe] is ∼1 dex higher for s18.8 than for b15-7. In both cases, the amount of ⁶⁴Zn produced from proton-rich conditions is comparable. They differ in the amount of more neutron-rich Zn isotopes. However, the abundances of A > 80 nuclei are somewhat uncertain as our tracer resolution is too coarse to adequately resolve the late neutrino-driven ejecta. For intermediate and light elements, any differences in the yields are mostly due to differences in the progenitor properties. For example, b15-7 ejects less ²⁸Si and other nuclei with $20\leqslant A\leqslant 40$ than s18.8. In return, the abundances of oxygen and carbon are somewhat higher in b15-7 when compared to s18.8. Other potential pre-explosion models of similar compactness as b15-7 and s18.8 are expected to yield similar explosion properties and nucleosynthesis yields for the iron group. Some differences in the lighter elements may be found, as these depend more strongly on the pre-collapse evolution than the iron group which is synthesized during the explosion.

In both cases (b15-7 and s18.8) our predicted yield of ⁴⁴Ti is too low when compared to observations, even when introducing a mixing-fallback prescription. To determine the extent of the mixing, multi-dimensional modeling is required [76–78]. In multi-dimensional simulations, the ejecta will experience convective overturn which can increase the amount of ⁴⁴Ti ejected. For example, Harris et al [79] found that the in situ yield of ⁴⁴Ti is greater than the yield obtained by postprocessing of two-dimensional CCSN simulations. A recent three-dimensional simulation found that ⁴⁴Ti can be synthesized in sufficient quantity and consistent with the observed spatial asymmetry even without invoking rapid rotation or a jet-driven explosion [80]. This indicates that the synthesis of ⁴⁴Ti could be a true multi-dimensional effect that cannot be captured in spherical symmetry.

A major open question is the compact remnant of SN 1987A. The observed neutrino signal from SN 1987A is consistent with the formation of a neutron star (NS), however it does not constrain the later evolution of the compact object (for example, continued accretion may lead to further collapse to a black hole). The detection of an electromagnetic signal from the compact object is hindered by the light from the ejecta and circumstellar medium [81]. To date, no compact object nor any evidence for a pulsar or other energy source have been detected despite searches across the electromagnetic spectrum. A recent study combined multi-wavelength observations and concluded that the most likely scenario is a dust-obscured neutron star that is thermally emitting [81]. We predict the mass of the neutron star born in the explosion for both models: b15-7 results in a cold neutron star of 1.48M_⊙ gravitational mass; s18.8 yields a cold neutron star of 1.41 M_⊙ gravitational mass. Hopefully, future observations will shed light on the nature and mass of the compact remnant of SN 1987A and test our predictions.

This work has been supported by the United States Department of Energy, Office of Science, Office of Nuclear Physics (award numbers SC0010263 and DE-FG02-02ER41216) and by the Research Corporation for Science Advancement. AH has been supported, in part, by a grant from the Science and Technology Commission of Shanghai Municipality (Grants No.16DZ2260200), the National Natural Science Foundation of China (Grants No.11655002), and the Australian Research Council through a Future Fellowship (FT120100363). GSI Darmstadt and the University of Basel are members in the European COST Action ChETEC (Chemical Elements as Tracers of the Evolution of the Cosmos, http://www.chetec.eu/).