# Nucleosynthesis and Clump Formation in a Core Collapse Supernova

## 1. Introduction

Overwhelming observational evidence (cf. references in Müller 1998) suggests that large-scale mixing processes took place in SN 1987 A and transported newly synthesized $`{}_{}{}^{56}\mathrm{Ni}`$ from its creation site near the collapsed core all the way out to the hydrogen envelope. Spectroscopic studies of SN 1987 F, SN 1988 A, SN 1993 J (Spyromilio, 1994, and references therein) and SN 1995 V (Fassia et al., 1998) indicate that such mixing is probably generic in core collapse supernovae. Indeed, artificial mixing of the radioactive ejecta within the (helium) envelope is indispensable in order to reproduce the light curves and spectra of Type Ib explosions using one-dimensional (1D) hydrodynamic models (Shigeyama, Nomoto, & Tsujimoto 1990; Woosley & Eastman 1997, and references therein).

These findings have instigated theoretical work on multidimensional supernova models focusing either on the role of convective instabilities in the delayed, neutrino-driven explosion mechanism within about the first second of evolution (Mezzacappa et al. 1998; Janka & Müller 1996; Burrows, Hayes, & Fryxell 1995; Herant et al. 1994; Miller, Wilson, & Mayle 1993), or on the growth of Rayleigh-Taylor instabilitites during the late evolutionary stages (Nagataki, Shimizu, & Sato 1998; Herant & Benz 1992; Müller, Fryxell, & Arnett 1991; Yamada & Sato 1991; Hachisu et al. 1990). However, multidimensional simulations which follow the evolution of the supernova shock from its revival due to neutrino heating, until its emergence from the stellar surface have not yet been performed.

In this Letter, we report on preliminary results of high-resolution two-dimensional (2D) computations which for the first time cover the neutrino-driven initiation of the explosion, the accompanying convection and nucleosynthesis as well as the Rayleigh-Taylor mixing within the first $`300`$ seconds of evolution.

## 2. Numerical method and initial data

We split our simulation into two stages. The early evolution ($`t1`$ s) is followed with a version of the HERAKLES code (T. Plewa & E. Müller, in preparation) which solves the multidimensional hydrodynamic equations using the direct Eulerian version of the Piecewise Parabolic Method (Colella & Woodward, 1984) augmented by the Consistent Multifluid Advection scheme of Plewa & Müller (1999) in order to guarantee exact conservation of nuclear species. We have added the input physics described in Janka & Müller (1996) (henceforth JM96) with the following modifications. General relativistic corrections are made to the gravitational potential following Van Riper (1979). A 14-isotope network is incorporated in order to compute the explosive nucleosynthesis. It includes the 13 $`\alpha `$-nuclei from $`{}_{}{}^{4}\mathrm{He}`$ to $`{}_{}{}^{56}\mathrm{Ni}`$ and a representative tracer nucleus which is used to monitor the distribution of the neutrino-heated, neutron-rich material and to replace the $`{}_{}{}^{56}\mathrm{Ni}`$ production when $`Y_\mathrm{e}`$ drops below $`0.49`$ (cf. Thielemann, Nomoto, & Hashimoto, 1996). Our initial data are taken from the 15 $`\mathrm{M}_{}`$ progenitor model of Woosley, Pinto, & Ensman (1988) which was collapsed by Bruenn (1993). The model is mapped to a 2D grid consisting of 400 radial zones ($`3.17\times 10^6\mathrm{cm}r1.7\times 10^9`$ cm), and 180 angular zones ($`0\theta \pi `$; cf. JM96 for details). A random initial seed perturbation is added to the velocity field with a modulus of $`10^3`$ of the (radial) velocity of the post-collapse model. The computations begin 20 ms after core bounce and are continued until 885 ms when the explosion energy has saturated at $`1.48\times 10^{51}`$ erg (this value has still to be corrected for the binding energy of the outer envelope). We will henceforth refer to this calculation as our “explosion model”.

The subsequent shock propagation through the stellar envelope and the growth of Rayleigh-Taylor instabilities is followed with the AMRA Adaptive Mesh Refinement (AMR) code (T. Plewa & E. Müller, in preparation). Neutrino physics and gravity are not included in the AMR calculations. Both do not influence the propagation of the shock during late evolutionary stages. However, gravity is important for determining the amount of fallback, a problem which is outside the scope of the present study. The equation of state takes into account contributions from photons, non-degenerate electrons, $`\mathrm{e}^+\mathrm{e}^{}`$-pairs, $`{}_{}{}^{1}\mathrm{H}`$, and the nuclei included in the reaction network. The AMR calculations are started with the inner and outer boundaries located at $`r_{\mathrm{in}}=10^8`$ cm (i.e. inside the hot bubble containing the neutrino-driven wind) and $`r_{\mathrm{out}}=2\times 10^{10}`$ cm, respectively. No further seed perturbations are added. Our maximum resolution is equivalent to that of a uniform grid of $`3072\times 768`$ zones. We do not include the entire star but expand the radial extent of the grid by a factor of 2 to 4 whenever the supernova shock is approaching the outer grid boundary, which is moved from its initial value out to $`r_{\mathrm{out}}=1.1\times 10^{12}`$ cm at $`t=300`$ s. Reflecting boundary conditions are used at $`\theta =0`$ and $`\theta =\pi `$ and free outflow is allowed across the inner and outer radial boundaries.

## 3. Results

The general features of our explosion model are comparable to the models of JM96. The most important difference is caused by our use of general relativistic corrections to the gravitational potential. Since the shock has to overcome a deeper potential well than in the Newtonian case, the luminosities (which are prescribed at the inner boundary) required to obtain a certain final explosion energy are roughly 20% higher than those of JM96. In particular, we have adopted the following set of parameters: $`L_{\nu _e}^0=2.8125\times 10^{52}\mathrm{erg}/\mathrm{s},L_{\nu _x}^0=2.375\times 10^{52}\mathrm{erg}/\mathrm{s},\mathrm{\Delta }Y_l=0.0875,\mathrm{\Delta }\epsilon =0.0625`$ (cf. JM96). The neutrino spectra and the temporal decay of the luminosity are the same as in JM96.

For the chosen neutrino luminosities the shock starts to move out of the iron core almost immediately. Convection between shock and gain radius sets in $`30`$ ms after the start of the simulation in form of rising blobs of heated, deleptonized material (with $`Y_\mathrm{e}0.5`$) separated by narrow downflows with $`Y_\mathrm{e}0.49`$. The shock reaches the Fe/Si interface at $`r=1.4\times 10^8`$ cm after $`100`$ ms. At this time the shocked material is still in nuclear statistical equilibrium and is composed mainly of $`\alpha `$-particles and nucleons. When the temperature right behind the shock drops below $`7\times 10^9`$ K, $`{}_{}{}^{56}\mathrm{Ni}`$ starts to form in a narrow shell. During the ongoing expansion and cooling, nickel is also synthesized in the convective region. However, this synthesis proceeds exclusively in the narrow downflows which separate the rising bubbles and have a sufficiently high electron fraction $`Y_\mathrm{e}`$. The convective anisotropies therefore lead to a highly inhomogeneous nickel distribution (Fig. 1).

Freezeout of complete silicon burning occurs at $`t250`$ ms, and convection ceases at $`t400`$ ms, at which time the flow pattern becomes frozen in. Subsequently, the entire post-shock region expands nearly uniformly. The post-shock temperature drops below $`2.8\times 10^9`$ K at $`t=495`$ ms, when the shock is about to cross the Si/O interface. Thus, our model shows only moderate oxygen burning (due to a non-vanishing oxygen abundance in the silicon shell), and negligible neon and carbon burning. This is caused by the specific structure of the progenitor model of Woosley et al. (1988) and may change when different (e.g. more massive) progenitors are used. In total, $`0.052\mathrm{M}_{}`$ of $`{}_{}{}^{56}\mathrm{Ni}`$ are produced, while $`0.10\mathrm{M}_{}`$ of material are synthesized at conditions with $`Y_\mathrm{e}<0.49`$ and end up as neutron-rich nuclei. Fig. 1 shows the $`{}_{}{}^{56}\mathrm{Ni}`$ distribution at $`t=885`$ ms. The density ratio between the dense regions which contain the nickel and the low-density, deleptonized material in the bubbles deeper inside is $`2.5`$.

During the next seconds the shock detaches from the formerly convective shell that carries the products of explosive nucleosynthesis and crosses the C+O/He-interface, initially accelerating and subsequently decelerating due to the varying density gradient. Twenty seconds after core bounce this unsteady propagation speed of the shock has led to a strong compression of the metal-containing shell. At the inner boundary of the shell pressure waves have steepened into a reverse shock. Rayleigh-Taylor instabilities start to grow at the Ni+Si/O- and C+O/He-interfaces and are fully developed and interact with each other at $`t=50`$ s. Nickel and silicon are dragged upward into the helium shell in rising mushrooms on angular scales from $`1^{}`$ to about $`5^{}`$, whereas helium is mixed inward in bubbles. Oxygen and carbon, located in intermediate layers of the progenitor, are swept outward as well as inward in rising and sinking flows. At $`t=300`$ s the densities and radial velocities between the dense mushrooms and clumps and the ambient medium differ by factors up to 5 and 1.3, respectively. As Fig. 2 shows, the fastest mushrooms have already propagated out to more than half the radius of the He core. The excessive outflows along the symmetry axis in Fig. 2 are a numerical artifact caused by the singularity of the polar axis in spherical coordinates. In order not to introduce large errors, only an angular wedge with $`15^{}\theta 165^{}`$ was used in all the analyses.

The remarkable efficiency of the mixing is clearly visible from Fig. 3 where the original onion-shell structure of the progenitor’s metal core has disappeared after 300 s. We observe large composition gradients between different regions of a mushroom and also between different mushrooms. Some of them contain $`{}_{}{}^{56}\mathrm{Ni}`$ mass fractions of more than 70% whereas others have only nickel admixtures of 20% or less but show high concentrations of silicon and oxygen. Composition contrasts of at least this magnitude have recently been observed in different filaments of the Cas A supernova remnant by Hughes et al. (2000). We find that also more than $`0.04\mathrm{M}_{}`$ of matter consisting of neutron-rich nuclei are mixed into the $`{}_{}{}^{56}\mathrm{Ni}`$ and $`{}_{}{}^{28}\mathrm{Si}`$ clumps.

The reverse shock substantially decelerates the bulk of the nickel from $`\mathrm{15\hspace{0.17em}000}`$ km/s at $`t=885`$ ms to 3200 – 4500 km/s at $`t=50`$ s, at which time the maximum velocities $`v_{\mathrm{Ni}}^{\mathrm{max}}`$ have dropped to 5800 km/s. After another 50 s, however, the clumps start to move essentially ballistically through the helium core and only a slight drop of $`v_{\mathrm{Ni}}^{\mathrm{max}}`$ from $`5000`$ km/s at $`t=100`$ s to $`4700`$ km/s at $`t=300`$ s occurs, when most of the $`{}_{}{}^{56}\mathrm{Ni}`$ has velocities below $`3000`$ km/s (Fig.4).

We have recently accomplished to follow the subsequent evolution of our model up to 16 000 s after core bounce. Similar to the situation at the C+O/He interface, the supernova shock leaves behind a dense (Rayleigh-Taylor unstable) shell at the He/H interface. While the entire shell is rapidly slowed down, a second reverse shock forms at its inner boundary (Fig. 2). Our high-resolution simulations reveal a potentially severe problem for the mixing of heavy elements into the hydrogen envelope of Type II supernovae like SN 1987A. We find that the fast nickel containing clumps, after having penetrated through this reverse shock, dissipate a large fraction of their kinetic energy in bow shocks created by their supersonic motion through the shell medium. This leads to their deceleration to $`2000`$ km/s in our calculations, a negative effect on the clump propagation which has not been discussed previously. A detailed analysis of this result will be subject of a forthcoming publication.

## 4. Conclusions

Our calculations provide a natural explanation for the mixing of $`{}_{}{}^{56}\mathrm{Ni}`$ into the helium layer, which is required to reproduce the light curves of Type Ib supernovae (Shigeyama et al., 1990). In this sense they justify the rather large seed perturbations ($`5\%`$) which were imposed by Hachisu et al. (1994) on the radial velocity given from spherically symmetric models of exploding helium cores at about 10 s after shock formation. In addition, they may have interesting implications for the modeling of Type Ib spectra (Woosley, private communication; compare Fig. 3 of this work with Fig. 7 in Woosley & Eastman 1997).

Our simulations suggest that ballistically moving, metal-rich clumps with velocities up to more than $`4000`$ km/s are ejected during the explosion of Type Ib (and Ic) supernovae. In case of Type II supernovae, however, the dense shell left behind by the shock passing the boundary between helium core and hydrogen envelope, causes a substantial deceleration of the clumps. The high iron velocities observed in SN 1987 A, for example, can therefore not be accounted for by our models. This problem would be reduced if the density profile at the He/H interface were smoother, leading to less strong variations of the shock velocity. Could this direct to a possible common envelope phase or merger history of the progenitor star (Hillebrandt & Meyer 1989; Podsiadlowski, Joss, & Rappaport 1990)? Alternatively, the high iron velocities in SN 1987 A could require additional energy input from Ni decay (Herant & Benz, 1992) or could imply a large global anisotropy of the explosion, e.g. associated with jets emerging from the collapse of a rapidly rotating stellar core (Fryer & Heger, 1999). Though final conclusions require three-dimensional calculations, our simulations indicate that all computations of Rayleigh-Taylor mixing in Type II supernovae carried out so far (including the case of SN 1987 A; Herant & Benz 1992) have been started from overly simplified initial conditions since they have neglected clump formation within the first minutes of the explosion.

The mixing of nucleosynthesis products as seen in our models poses another potential problem. A significant fraction ($`0.04\mathrm{M}_{}`$ of a total of $`0.1\mathrm{M}_{}`$) of the neutron-rich nuclei which are synthesized in regions with $`Y_\mathrm{e}<0.49`$ is dragged outward by the nickel containing mushrooms and clumps. At least this amount will most likely be ejected in the explosion. A detailed analysis of the nuclear composition is necessary to tell whether this is in conflict with limits from Galactic chemical evolution models, which allow for at most $`10^3`$$`10^2\mathrm{M}_{}`$ of neutron-rich nuclei to be ejected per supernova event (e.g., Herant et al. 1994; Thielemann et al. 1996). Later fallback will not solve this potential problem: How could it disentangle the clumpy nickel ejecta from their undesirable pollution? A better knowledge of the luminosities and spectra of electron neutrinos and antineutrinos emitted from the nascent neutron star is therefore needed to perform more reliable calculations of the neutronization of the neutrino-heated ejecta.

During the computations we became aware of oscillations with angle in parts of the postshock flow (Figs. 1 and 2). These are caused by the “odd-even decoupling” phenomenon associated with grid-aligned shocks (Quirk, 1994). As a consequence, the maximum nickel velocities, $`v_{\mathrm{Ni}}^{\mathrm{max}}`$, obtained in our AMR calculations have probably been overestimated by $`25\%`$ because the growth of some of the mushrooms was influenced by the perturbations induced by this numerical defect. The main results of our study, however, are not affected. We note that a large number of supernova calculations performed with codes based either on the direct Eulerian (cf. Figs. 22 and 24 in Burrows et al. 1995, Fig. 20 in JM96), or the Lagrangean with remap formulation of the PPM scheme (Mezzacappa et al., 1998) seem to be affected by this numerical flaw. We defer a detailed analysis of this problem to a forthcoming publication.

We thank S. Woosley for profiles of the progenitor star and S. Bruenn for his post-bounce core model. We acknowledge support by P. Cieciela̧g and R. Walder concerning visualization and by the crew of the Rechenzentrum Garching where the simulations were performed on the NEC SX-4B and CRAY J916 computers. TP was supported by grant 2.P03D.004.13 from the Polish Committee for Scientific Research, HTJ by DFG grant SFB-375 für Astro-Teilchenphysik.