THREE-DIMENSIONAL SIMULATIONS OF MIXING INSTABILITIES IN SUPERNOVA EXPLOSIONS

We use a computational grid in spherical polar coordinates consisting of 1200 × 180 × 360 zones in r-, θ-, and ϕ-direction. The equidistant angular grid has a resolution of 0935 in θ-direction and 1° in ϕ-direction, covering the whole sphere except for a (double) cone with a half opening angle of 58 around the symmetry axis of the coordinate system (i.e., 0.0325 π ⩽ θ ⩽ 0.9675 π, and 0 ⩽ ϕ ⩽ 2π). Reflective boundary conditions are imposed in θ-direction, and periodic ones in ϕ-direction. The clipping of the computational grid around the symmetry axis avoids a too restrictive CFL time step condition. We have not observed any numerical artifacts as consequences of this grid constraint.

The radial mesh is logarithmically spaced between a time-dependent inner boundary, initially located at a radius of 200 km, and a fixed outer boundary at 3.9 × 10¹² cm. It has a maximum resolution of 2 km at the inner boundary, and a resolution of 4 × 10¹⁰ cm at the outer one (corresponding to a roughly constant relative radial resolution Δr/r ≈ 1%). We allow for free outflow at the outer boundary, and impose a reflective boundary condition at the inner edge of the radial grid.

During the simulation we move (approximately every 100th time step) the inner radial boundary to larger and larger radii to relax the CFL condition. This cutting of the computational grid reduces the number of (the logarithmically spaced) radial zones from 1200 at the beginning to about 400 toward the end of the simulation, but involves only the innermost few percent of the initial envelope mass. We convinced ourselves by means of 2D test calculations that this removal of mass has no effect on the dynamics and mixing occurring at larger radii (the cut radius exceeds at no time a few percent of the radius of the Rayleigh–Taylor finger tips).

The initial model used for our 3D simulations is based on 3D hydrodynamic explosion models of Scheck (2007), who followed the onset and early development of neutrino-driven explosions in 3D with the same code, numerical setup, and input physics as described in all details for 1D and 2D simulations in Scheck et al. (2006). These runs were started from 1D stellar collapse models a few milliseconds after core bounce. The explosion was initiated and powered by the energy input of neutrinos in the post-shock layer, where vigorous convective activity was present and supported the revival of the stalled prompt shock. A defined value of the explosion energy was obtained by suitably choosing the magnitude of the neutrino luminosities that were imposed at the retreating inner grid boundary, which replaced the excised, contracting dense core of the forming neutron star.

The stellar progenitor considered here is a spherically symmetric model of a 15.5 M_☉ blue supergiant with a radius of 4 × 10¹² cm (Woosley et al. 1988, S. E. Woosley 1999, private communication; see also Kifonidis et al. 2003). Since the models of Scheck (2007) included only the central part of the star out to the middle of the C/O shell at r = 1.7 × 10⁹ cm, the envelope of the progenitor had to be added for the long-time simulations presented here.

The 15.5 M_☉ blue supergiant stellar progenitor model was evolved by Woosley et al. (1988) using a nuclear reaction network including a large set of nuclei. However, Scheck (2007) did not consider nuclear burning, and he included only a small number of nuclear species (p, n, α, and $^{54}\rm Mn$ representing the heavy element fraction) in their simulations to save computational costs. For our initial model we kept the nuclear composition of the 3D explosion model of Scheck (2007) inside the shock radius, redefining $^{54}\rm Mn$ into $^{56}\rm Ni$. The nuclear composition at radii greater than the shock radius was taken from the 15.5 M_☉ progenitor model.

Scheck (2007) simulated one 3D hydrodynamic explosion model with 2° angular resolution, which reached a final time of 0.58 s after bounce, and another one with 3° angular resolution, which covered an evolution time of 1 s after core bounce. Although both models were exploited for our study, we report here only on simulations using the former model that according to the data of Scheck (2007) has an (unsaturated) explosion energy of 0.6 × 10⁵¹ erg. The deformation of the shock front at t ∼ 0.5 s after core bounce and the asymmetry of the entropy distribution behind the shock in the form of rising, inflating hot plumes and infalling, narrow, cool downflows are shown in Figure 1. The structure of the explosion at this time looks very similar to comparable results by Fryer & Warren (2002). Besides following the shock's propagation through the stellar envelope of this initial model with our 3D code, we also performed a 3D simulation where we artificially increased the explosion energy of this initial model to 1.0 × 10⁵¹ erg by extending the thermal energy input as accomplishable by ongoing neutrino heating.

Zoom In Zoom Out Reset image size

We used a tabulated equation of state (EOS; Timmes & Swesty 2000) that considers contributions of an arbitrarily degenerate and relativistic electron–positron gas, of a photon gas, and of a set of ideal Boltzmann gases, consisting of the eight nuclear species (n, p, α, $^{12}\rm C$, $^{16}\rm O$, $^{20}\rm Ne$, $^{24}\rm Mg$, and $^{56}\rm Ni$) included in our initial model. Note that the pre-shock layers contained only a small region of silicon at the time the explosion model of Scheck (2007) was adopted as initial data for our long-time SN runs. Since nuclear burning was not taken into account in our simulations, but it is possible that this silicon will be explosively burned to nickel, we lumped the small amount of silicon and the iron-group elements in the post-shock layer together to what we followed as $^{56}\rm Ni$ in our hydrodynamic calculations.

We performed a set of thirteen 2D and three 3D simulations and compared the results obtained in three dimensions to those of the corresponding 2D models. In this paper, we focus our discussion on one 3D model that has an explosion energy of 10⁵¹ erg. The main findings for this particular 3D model are qualitatively similar to those we obtained for less energetic 3D explosions, but the resulting maximum values for the asymptotic clump velocities vary with the explosion energy roughly as $v_{\mathrm{clump}}\propto \sqrt{E_\mathrm{exp}}$.

The 2D simulations reported here were started from different meridional slices of the 3D initial explosion model. This ensures as much similarity of the initial conditions as possible and the same numerical resolution for the 2D and 3D runs. But the differences between the chosen slices give rise to some variation of the initial asymmetry and explosion energy among the various 2D models, which manifests itself both in different clump structures and velocities, and the amount of mixing in the stellar envelope.

A set of 13 meridional 2D slices was placed at azimuthal angles of 4°, 36°, 76°, 108°, 128°, 148°, 180°, 220°, 228°, 252°, 292°, 324°, and 360°. Three of these slices (at 128°, 220°, and 228°) were chosen at locations where the growth of Rayleigh–Taylor fingers was observed to be particularly strong in 3D. This guaranteed that such regions were not missed by our comparison of 2D and 3D calculations. It turned out that the largest 3D Rayleigh–Taylor structures at late stages develop in the directions of the biggest initial buoyant bubbles, expanding toward the bottom edge and upper right corner (at six o'clock and one o'clock, respectively) of Figure 1.

Both in the 3D and 2D simulations we neglected the influence of gravity on the motion of the ejecta. While this has no important impact on the dynamics of the expanding ejecta, the amount of fallback is underestimated. However, that way we could avoid the accumulation of mass near the inner (reflecting) radial boundary which would have implied a considerably more restrictive CFL condition.

The dynamical evolution of the explosion after its launch by neutrino energy deposition and the growth of instabilities in the considered 15.5 M_☉ progenitor were described in much detail by Kifonidis et al. (2003). It was shown there that the asymmetries associated with the convective activity developing in the post-shock region during the neutrino-heating phase act as seeds of secondary RTIs at the composition interfaces of the exploding star. At about 100 s dense Rayleigh–Taylor fingers containing the metals (C, O, Si, iron-group elements) have grown out of a compressed shell of matter left behind by the shock passing through the Si/O and (C+O)/He interfaces. These fingers grow quickly in length and while extending into the helium shell, fragment into ballistically moving clumps and filaments that propagate faster than the expansion of their environment.

While for a 2D model with explosion energy around 1.8×10⁵¹ erg the Si and Ni containing structures still move with nearly 4000 km s⁻¹ (oxygen has velocities up to even 5000 km s⁻¹) at 300 s, a strong deceleration occurs when the metal clumps enter the relatively dense He shell that forms after the shock passage through the He/H interface. At t ≳ 10,000 s the metal carrying clumps have dissipated their excess kinetic energy and propagate with the same speed as the helium material in their surroundings. In the presence of a hydrogen envelope and the corresponding deceleration of the shock as it propagates through the inner regions of the hydrogen layer (in which case the helium "wall" below the He/H interface builds up), Kifonidis et al. (2003) could not observe any metal clumps that achieve to penetrate into the hydrogen envelope, in contrast to what was observed in the case of SN 1987A.

The 2D calculations performed in the course of the present work confirm these findings (when the velocities are appropriately rescaled with $\sqrt{E_\mathrm{exp}}$ to account for the lower explosion energies of E_exp ≲ 10⁵¹ erg considered here compared to roughly twice this value employed by Kifonidis et al. 2003). The evolution and the final results for the mass distributions of different chemical elements in velocity space and mass space are in good quantitative agreement with the models of Kifonidis et al. (2003), although the latter had much superior numerical resolution because of the use of a sophisticated adaptive mesh refinement algorithm. We are therefore confident that at least with respect to gross (quasi-1D) features, the results presented here are numerically converged, even if many of the structural details, including plumes and bullets, may not be numerically converged.

Kifonidis et al. (2006), investigating explosions also with E_exp ∼ 2 × 10⁵¹ erg, proposed a cure of the metal-H mixing problem by invoking a sufficiently large dipolar or quadrupolar asymmetry of the beginning explosion with a globally aspherical shock wave. This was motivated by recent hydrodynamical studies that show that large-amplitude, low-ℓ mode oscillations (ℓ defining the order of an expansion in spherical harmonics) due to the SASI can be an important ingredient in the explosion mechanism and in the sequence of processes that leads to observed asymmetries of SN explosions and neutron star kicks (see, e.g., Blondin et al. 2003; Blondin & Mezzacappa 2007; Scheck et al. 2004, 2006; Scheck 2007; Scheck et al. 2008; Burrows et al. 2006, 2007; Marek & Janka 2009). This on the one hand led to higher initial maximum velocities of the metal-rich clumps and thus their faster propagation through the He core on a timescale shorter than the build-up of the thick He "wall" that prevented their penetration into the hydrogen shell in the older calculations. On the other hand it triggered the growth of RMI at the He/H interface due to the deposition of vorticity by the shock. This led to strong mixing of hydrogen inward in mass space and down to low velocities, an effect that has turned out to be necessary to explain the shape of the light curve maximum in the case of SN 1987A (Utrobin 2004).

These conclusions, however, were based on 2D models. Therefore the question remained whether 3D might yield differences and whether also in 3D one would have to advocate a pronounced global low-ℓ mode asphericity of the beginning explosion as an explanation of the high velocities (≳3000 km s⁻¹) of nickel clumps in the hydrogen shell and of the strong inward mixing of hydrogen observed in SN 1987A. The 3D simulations performed in the present work demonstrate that this is not the case and the mentioned observational features can be accounted for by 3D models even with explosion energies around 10⁵¹ erg for the considered progenitor star.

In order to quantify the global shock deformation of our 3D initial models compared to the 2D cases investigated by Kifonidis et al. (2006), we did not perform an analysis in terms of an expansion in spherical harmonics modes. Instead, sizeable differences become obvious already when a tri-axial ellipsoid is constructed as envelope of the bumpy shock surface obtained by Scheck (2007) such that the major axis of the ellipsoid coincides with the direction of the strongest shock expansion. While the 3D model used for the calculations of the present paper has an axis ratio of 1.04:1.02:1.00 at t = 0.58 s after bounce—another 3D explosion run of Scheck (2007) yields 1.16:1.06:1.00 at t = 1 s—the shocks of the 2D configurations studied by Kifonidis et al. (2006) typically had much larger axis ratios of about 1.5:1.0.

The deformed shock causes strong RMI to develop at the (C+O)/He and He/H interfaces in the 2D models discussed here. Compared to the runs of Kifonidis et al. (2006), however, the smaller asphericity of the shock leads to a much slower growth of the instability. In our 2D runs the corresponding vortices at the He/H boundary at 13,000 s have reached a size relative to the average interface radius similar to what was observed at 3000 s in the models of Kifonidis et al. (2006). We note in passing that apart from resolution-dependent fine structure, the results of the latter paper could be very well reproduced concerning the RMI vortex shape, size, and locations by simulations performed for the same initial conditions but with the grid setup and resolution used in the present work (see Hammer 2009). The growth of the RMI observed in average cases of our 2D calculations with only weak low-ℓ mode shock deformation—in agreement with analytic growth rate estimates; Hammer (2009)— is much too slow to produce any significant extent of hydrogen–helium mixing so that the final outcome of our 2D models basically confirms the findings of Kifonidis et al. (2003).¹ In the 3D models RMI distortions can be seen at the (C+O)/He interface and are likely to contribute to the turbulent mixing of the metal core with the helium shell of the exploding star. At the H/He interface, however, where the shock is very close to being spherical, no clear RMI activity becomes visible before it is penetrated by fast, metal-carrying clumps that have been able to pass through the helium layer with still high velocities.

Figure 2 displays the development of these fast-moving clumps during our 3D explosion run by showing surfaces of constant mass fractions of carbon, oxygen, and nickel for two different viewing directions and two different times (350 s and ∼9000 s after core bounce). Figure 3 provides a volume-rendered image of the composition distribution at the later time, while Figure 4 gives composition information on cut planes through the mixed stellar core and some of the major plumes of different types. Finally, Figures 6 and 7 present normalized mass distributions of various nuclear species in the radial-velocity and enclosed-mass space for our 3D simulation compared to an "average" 2D result at several times after core bounce, and Figure 8 provides information for the spread of the hydrogen and nickel mass distributions in our set of 2D runs at the end of the simulations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We stress that what we denote as "nickel" here and in the following actually includes the contributions of silicon and of the iron-group elements from the post-shock region. This makes sense because the shock front in the initial model has nearly reached the oxygen layer (which implies that only a very small fraction of the original silicon is left) and because we do not account for nuclear burning in our simulations (for which reason no silicon is produced later on); see Section 2.3. In the 3D model the total mass of free neutrons is 6 × 10⁻⁴ M_☉, of hydrogen 8.16 M_☉, of helium 5.40 M_☉, of carbon 0.12 M_☉, of oxygen 0.20 M_☉, of neon 0.048 M_☉, of magnesium 0.005 M_☉, and of "nickel" 0.212 M_☉. About 1.35 M_☉ of the 15.5 M_☉ progenitor are thus absorbed into the neutron star. In the 2D runs based on selected meridional slices, the ejecta masses of the different nuclear species can differ insignificantly from the listed numbers.

Figure 2 shows that at 350 s oxygen- and nickel-rich, dense fingers, some of which display the typical mushroom shape of Rayleigh–Taylor structures, penetrate through the outer boundary of the carbon layer of the exploding star. The corresponding ¹²C surface is defined by a mass fraction of 3% (note that also exterior to this surface, in the surrounding helium and hydrogen layers, a small but non-zero mass fraction of carbon is still present and representative of the heavy elements that define the metallicity of the progenitor star). The most prominent features at that time have been seeded by the biggest high-entropy bubbles visible in Figure 1 and develop to the largest mushrooms visible more than 8000 s later.

The chemical composition of the various kinds of extended plumes in this late stage is very interesting. All of them contain high mass fractions of hydrogen and helium, which account for about 50%–70% of their mass with spatial variations in individual filaments. The dominant heavy species, however, differs between different structures. The largest mushrooms carry predominantly nickel and may contain a narrow spine with oxygen. In contrast, many of the far more numerous, less extended and more narrow fingers are composed of a mixture of nickel and oxygen. Also clumps with a dominant mass fraction of oxygen can be identified. More useful for such an analysis than the surface images of Figure 2 is the visualization by volume-rendering techniques in Figure 3. While the surfaces of constant mass fraction can be misleading, because their shape depends strongly on the chosen value and no information about the interior composition is provided, the ray-tracing image shows that the largest Rayleigh–Taylor mushrooms carry mostly nickel, the majority of the more extended fingers are composed of a mix of elements, and only some clumpy, relatively compact, knotty features (clearly visible along the edge of the pastel-colored and whitish glowing core in Figure 3) contain a dominant fraction of oxygen. This is confirmed by the planar cuts through two typical Rayleigh–Taylor structures in Figure 4. The upper left panel of this figure, in which a 2D cut through the whole core is displayed, demonstrates that on the global scale there are extended and separated regions where different species are the most abundant chemical constituent. In the core the hydrogen and helium admixture is much smaller than in the extended fingers and plumes. Except in some smaller "bubbles," where H+He dominate, heavies contribute 50%–80%, in some clumps even more than 90%, of the core composition.

Zoom In Zoom Out Reset image size

The time sequences of the normalized mass distributions of the nuclear species in the radial-velocity and enclosed-mass space² of Figures 6 and 7 visualize the gradual acceleration of the outer He and H layers of the exploding star by the outgoing shock front. Simultaneously, the initially fast metals of the core penetrate into the outer composition layers and are decelerated. Both the deceleration and the mixing depend crucially on the hydrodynamic instabilities that develop at the composition interfaces, seeded by the initial explosion asymmetries. At 350 s the helium "wall" can be discerned as the peak of the distribution between 3500 km s⁻¹ and 4500 km s⁻¹ at an enclosed mass between ∼2.5 M_☉ and ∼4.5 M_☉. From a comparison of the plots for 2D and 3D it is obvious that nickel and oxygen experience a much stronger deceleration during the subsequent evolution and far less efficient mixing into the hydrogen envelope in the 2D case (the plane of the meridional cut considered in Figures 6 and 7 is displayed in Figure 5).³ Since neon and magnesium are carbon-burning ashes and well mixed with oxygen in the oxygen layer of the progenitor star, the distributions of these three elements closely match each other during the mixing evolution and can be hardly distinguished on plots. We therefore have combined them into one line in Figures 6 and 7.

Zoom In Zoom Out Reset image size

A closer inspection of Figures 6 and 7 reveals that the metal distributions in both 2D and 3D models peak at the location of the metal core (enclosed mass M_r ≲ 2.5 M_☉) with similar velocities (around 2000–3000 km s⁻¹ at 350 s and 1000 s, around 1500–2500 km s⁻¹ at 2600 s, and around 1000 km s⁻¹ at 9000 s). In 3D, however, the nickel and oxygen distributions possess much more massive and wider high-velocity tails up to velocities of more than 6000 km s⁻¹ at 350 s. In the long run this is causal for the large discrepancies of the radial composition mixing. Two effects appear to contribute to the formation and continued presence of these high-velocity tails.

• 1.  
  In 3D the RTIs at the composition interfaces (in particular the (C+O)/He boundary) grow more rapidly. This leads to the formation of extremely fast (mostly nickel and less oxygen) clumps with velocities in excess of 5000 km s⁻¹ at 350 s (upper right panel of Figure 6). These dense bullets penetrate deep into the hydrogen layer even before the helium wall has formed (see the upper two right panels in Figure 7). The more rapid growth of the RTIs also seems to be the reason for the efficient mixing of hydrogen into the metal core, which is clearly different from the 2D case already at ∼1000 s.
• 2.  
  A larger fraction of the nickel and oxygen with a velocity of ∼4000 km s⁻¹ at 350 s, which is in the ballpark of the maximum nickel and oxygen velocities of the 2D models, experiences less velocity reduction in the 3D simulation. This is obviously linked to a different propagation and deceleration of the dense metal-containing clumps as they move quasi-ballistically through the helium shell and inner part of the hydrogen envelope under the influence of drag forces exerted by the surrounding medium. While the fastest bullets with a speed of 5000–6000 km s⁻¹ at 350 s still move with 3000–4500 km s⁻¹ at t>2600 s, the clumps with initially ∼4000 km s⁻¹ are decelerated to 2000–3000 km s⁻¹. In contrast, in 2D none of the material of the metal core has retained velocities of more than 2000 km s⁻¹ at the time when the shock reaches the stellar surface (t ∼ 8000 s). We will return to a detailed discussion of this very interesting phenomenon in Section 4, where we will present simple analytic arguments for an explanation.

Zoom In Zoom Out Reset image size

In Figure 8, the normalized mass distributions of hydrogen and nickel in the radial velocity space are plotted at t ∼ 9000 s shortly after the shock breakout from the stellar surface (when the expansion is close to become homologous) for our 3D run in comparison to the compilation of results obtained for all 2D simulations with different chosen meridional slices. While the maximum hydrogen velocities are very similar in 2D and 3D and reach up to about 6000 km s⁻¹, a clearly bigger hydrogen mass attains the lowest velocities (below 1000 km s⁻¹) in the 3D case. Correspondingly, a small fraction of the hydrogen (some hundredths of a solar mass) is mixed deep into the metal core of the star to regions with enclosed masses of M(r) ≲ 2 M_☉ (see also the right panels of Figure 7). These findings are very similar to the results obtained in 2D by Kifonidis et al. (2006; cf. their Figures 6, 8, and 9), who considered cases with largely aspherical beginning of the explosion as discussed in Section 3.1). In 3D similarly strong inward mixing of hydrogen does not require a very pronounced global asymmetry of the early blast. RMI at the shock-distorted He/H composition interface plays no crucial role in our 3D simulations, because the outgoing SN shock is significantly less deformed than in the 2D models discussed by Kifonidis et al. (2006), see Section 3.1.

Zoom In Zoom Out Reset image size

Even more impressive than the hydrogen differences is the large discrepancy of 2D and 3D with respect to the nickel distribution in velocity space and the maximum nickel velocities. While in 2D even the fastest nickel clumps move with just 2200 km s⁻¹ at t ∼ 9000 s (Figure 8), a large fraction of the nickel retains velocities of more than 3000 km s⁻¹ and some even more than 4000 km s⁻¹ in 3D (Figure 8). From Figure 6 we can conclude that similarly large 2D versus 3D differences are obtained for the oxygen distribution in velocity space, whereas the hydrogen and helium distributions of the 2D case are close to the results of the 3D run (except for the low-velocity tail of the hydrogen distribution mentioned above), and the carbon distribution exhibits slight differences only between about 1000 km s⁻¹ and 2500 km s⁻¹. Interestingly, a significantly larger fraction of the nickel than of the oxygen propagates with the highest speeds of clumps. This is compatible with the fact that the metal content of the biggest and most extended Rayleigh–Taylor structures in Figure 2 is dominated by nickel, while oxygen is concentrated mostly in thin, shorter cores of the mushrooms (see Figure 4).

It should also be noted that the spread of the maximum nickel velocities obtained in the set of 2D runs in Figure 8 is considerably smaller than the corresponding spread of the 2D results in the mass space. The maximum Ni velocities vary only between about 1700 km s⁻¹ for average and 2200 km s⁻¹ for extreme 2D cases (in which the meridional slices were placed at the locations of the biggest nickel mushrooms in 3D). In the former models nickel is transported to slightly more than 5 M_☉, whereas in the latter ones some of the nickel gets mixed out to an enclosed mass of even 8–9 M_☉. In all the cases, however, the mixing of iron-group elements into the hydrogen layer is much less strong in 2D than in 3D.

Zoom In Zoom Out Reset image size

Figure 7, showing an average 2D run, confirms these findings. In the latest stages nickel and oxygen can be observed at an enclosed mass of about 5 M_☉ (in rough agreement with the results of Kifonidis et al. 2003, see their Figure 15), which is close to the base of the hydrogen envelope.⁴ In typical 2D cases, the metal-containing clumps thus get stuck in the dense helium wall as pointed out by Kifonidis et al. (2003). In contrast, in the 3D run nickel and oxygen are mixed deep into the hydrogen layer (and also a small amount of carbon is carried from the metal core into the hydrogen shell). The fastest nickel and oxygen clumps are able to penetrate nearly to the surface of the exploding star (in the mass space).

The two humps that can be discriminated in the mass distribution of these nuclear species well ahead of the bulk of the metal core in the mass and velocity space at 9000 s, are linked to the few, but big and far-reaching structures, and to the greater number of smaller fingers that extend more uniformly to all directions in Figure 2. The bigger bullets contain about three to four times more nickel than oxygen. They carry in total a mass of ∼2×10⁻² M_☉ with velocities of 3000–4500 km s⁻¹ for the nickel and 3000–4000 km s⁻¹ for the oxygen. In contrast, the smaller finger-like structures move with 1500–3000 km s⁻¹ and contain a total mass of ∼4×10⁻² M_☉.

It is amazing that in our 3D simulations the fastest of the metal clumps have achieved to overtake most of the hydrogen at the time the shock breaks out (at roughly 8000 s). The nickel and oxygen mixing into the hydrogen envelope is clearly more efficient than even in the globally asymmetric explosions studied by Kifonidis et al. (2006), where at 10,000 s the metals were seen to be distributed only up to an enclosed mass of about 10 M_☉ (see their Figure 6), although the peak metal velocities were roughly the same at ∼300 s.

We refrain here from drawing far reaching conclusions on the observational consequences of this finding, e.g., concerning the early visibility of X-ray and γ-ray signatures of such strong mixing. While we consider the 2D/3D differences for the same explosion model and the same and fixed numerical resolution of our first explorative study as enlightening, we think that more simulations with different explosion energies, different progenitor stars, and in particular also with finer grid zoning are desirable before the results can be interpreted quantitatively with respect to observations of SNe.

A first qualitative understanding can be obtained by considering objects moving in a stellar medium with a time-independent density ρ(r). We assume that the clump mass m = ρ_cV_c remains constant and the volume V_c and density ρ_c of the clumps adjust to the local environmental conditions accordingly. With a velocity v relative to the surrounding gas the equation of motion of the clumps at a radius r can be written as

$$\begin{equation} m {{{d}}v(r)\over {{d}}t}= {1\over 2} \, m {{{d}}v^2(r)\over {{d}}r} = -F_{\mathrm{drag}}(r)\,, \end{equation} \tag{ 1 }$$

where the drag force for high Reynolds numbers (Re = ρvV_c(A_cμ)⁻¹ ≫ 1 with μ being the dynamical viscosity and A_c the area the object presents to the flow) is given by

$$\begin{equation} F_{\mathrm{drag}}(r) = {1\over 2}\, C A_\mathrm{c}(r) \rho (r) v^2(r) \,. \end{equation} \tag{ 2 }$$

The drag coefficient C is a dimensionless parameter of order unity, which depends on the shape and the surface structure of the moving object. Equation (1) with Equation (2) has an analytic solution for the velocity as a function of distance:

$$\begin{equation} {v^2(r)\over v^2_0}= \exp \left(-\,{C\over m}\int _{r_0}^r {{d}}r^{\prime }\, A_\mathrm{c}(r^{\prime }) \rho (r^{\prime }) \right)\,\equiv \,{e}^{-2 f(r)}\,, \end{equation} \tag{ 3 }$$

when v₀ is the initial velocity of the object at radial position r₀. The exponential decay of the initial velocity due to the influence of the drag force depends on the density of the medium through which the object propagates, on the drag coefficient, the mass of the object, and the area it presents to the flow.

In the following, we will discuss how the term in the exponent depends on the dimensionality of the simulation. Because of the assumed axial symmetry, in 2D all objects have a toroidal form. Let us consider in 2D an ideal torus with radius r_t, whose center has a distance r from the origin of the stellar grid and a distance rsin θ from the symmetry axis, and compare its outward propagation in radial direction with that of a sphere with radius r_s in 3D also at radial distance r from the center.

These objects, torus and sphere, are assumed to have a fixed shape but their size can vary as they move in radial direction from the initial position r₀ to r. In the beginning, the torus radius is r_t,0, and the spherical blob has a radius r_s,0. In order to compare the propagation of structures in 3D and 2D as in the hydrodynamical simulations, we suppose that the torus and sphere have initially the same cross section and thus radius: r_s,0 = r_t,0. Since we assume the same density inside the objects, ρ_t,0 = ρ_s,0, but the torus volume V_t,0 = 2π²r²_t,0r₀sin θ in general differs from the volume of the sphere, $V_{{\mathrm s},0} = {4\pi \over 3} r_{{\mathrm s},0}^3$, the masses of torus and sphere are also unequal: m_t ≠ m_s. Besides these different masses, the two objects also present different areas to the flow, A_t,0 = 4πr_t,0r₀sin θ and A_s,0 = πr²_s,0.

Moving in radial direction, the clumps are assumed to maintain pressure equilibrium with the surrounding medium. As a consequence, their density and therefore their volume and cross section will change. For the torus one has the scaling relations V_t = V_t,0(r_t/r_t,0)²(r/r₀) and A_t = A_t,0(r_t/r_t,0)(r/r₀), while for the sphere they are V_s = V_s,0(r_s/r_s,0)³ and A_s = A_s,0(r_s/r_s,0)². Using this and the equality of the initial density and radii of torus and sphere, the arguments in the exponential function of Equation (3) can be compared for the two geometries to yield the ratio:

$$\begin{eqnarray} {f_{\mathrm t}(r)\over f_{\mathrm s}(r)} &=& {C_{\mathrm t}\over C_{\mathrm s}}{8\over 3\pi } \left[ \int _{r_0}^r \! {d}r^{\prime }\,\rho (r^{\prime }) \,{r_{\mathrm t} \over r_{{\mathrm t},0}}\, {r^{\prime }\over r_0} \right] \nonumber\\ && \times \left[ \int _{r_0}^r \! {d}r^{\prime }\,\rho (r^{\prime }) \,\left({r_{\mathrm s} \over r_{{\mathrm s},0}} \right)^{\! 2} \right]^{-1}. \end{eqnarray} \tag{ 4 }$$

In this equation the radii of the torus, r_t, and of the sphere, r_s, in general depend on the distance r from the grid center, if the clump propagates through a stellar medium with a pressure gradient. It should be noted that the ratio f_t(r)/f_s(r) is nearly independent of both the mass difference and the initial difference of the torus and sphere areas perpendicular to the radial direction (the numerical factor 8/(3π) is very close to unity). Instead, the deceleration of the motion by drag forces depends crucially on the evolution of the areas as reflected by the radius-dependent functions in the integrands.

We assume that the pressure in the clump (torus or sphere) scales with the density in the clump as P_c ∝ ρ^m_c ∝ V^−m_c, and that the pressure in the stellar environment follows a radial profile P(r) ∝ r⁻ⁿ. Using this and the requirement of pressure equilibrium, we find that the radii of torus and sphere change with r according to

$$\begin{equation} {r_{\mathrm t}\over r_{{\mathrm t},0}} = \left({r\over r_0}\right)^{\! {1\over 2}{n\over m} - {1\over 2}} \,,\quad \ \ {r_{\mathrm s}\over r_{{\mathrm s},0}} = \left({r\over r_0}\right)^{\! n\over 3m}. \end{equation} \tag{ 5 }$$

Plugging Equation (5) into Equation (4), we obtain

$$\begin{eqnarray} {f_{\mathrm t}(r)\over f_{\mathrm s}(r)} &=& {C_{\mathrm t}\over C_{\mathrm s}}{8\over 3\pi } \left[ \int _{r_0}^r \! {d}r^{\prime }\,\rho (r^{\prime }) \, \left({r^{\prime }\over r_0} \right)^{\! {1\over 2}{n\over m} + {1\over 2}} \right] \nonumber\\ && \times \left[ \int _{r_0}^r \! {d}r^{\prime }\,\rho (r^{\prime }) \,\left({r^{\prime }\over r_0} \right)^{\! 2n\over 3m} \right]^{-1}. \end{eqnarray} \tag{ 6 }$$

The torus area perpendicular to the radial direction grows faster than that of the sphere if n < 3m, in which case the integral in the numerator increases more rapidly with r than the one in the denominator. This difference is a consequence of the fact that any radial motion leads to a larger distance of the torus from the polar axis of the grid. This stretches the torus and causes an adjustment of its cross section according to pressure equilibrium. Since in addition to the area effect also the drag coefficient of the toroidal object is likely to be larger than that of the sphere, C_t ≳ C_s, the motion of the torus is clearly more strongly damped by the drag than that of the sphere if n < 3m.

In the SN the clumps propagate through the medium of the dense shell of shock-accelerated progenitor matter. In this shell the density and pressure gradients are nearly flat, i.e., ρ ∼ const and n ∼ 0 (see, e.g., Nomoto et al. 1994, Figures 31 and 32). This case is particularly suitable to illustrate the inadequacy of axisymmetric toroidal clumps in 2D instead of spherical clumps in 3D. While the latter propagate through the homogeneous medium without any change of their area perpendicular to the direction of the motion, the tori are unavoidably stretched as they reach larger radii r, which leads to an increase of the cross-sectional area they present to the flow. Equation (6) yields

$$\begin{equation} {f_{\mathrm t}(r)\over f_{\mathrm s}(r)} \left({C_{\mathrm t}\over C_{\mathrm s}}{8\over 3\pi } \right)^{\! -1} = {2\over 3}\, \left[ \left({r\over r_0}\right)^{\! {3\over 2}} - 1 \right]\, \left[ {r\over r_0} -1 \right]^{-1}, \end{equation} \tag{ 7 }$$

which is larger than unity for r>r₀. A stronger deceleration of the toroidal clump is thus caused by the artificial constraints associated with the axisymmetry of all structures in the 2D geometry.

In a second step we intend to obtain a more quantitative insight by considering the clump propagation in the dense helium layer that forms after the passage of the outgoing shock through the He/H interface (for details, see Kifonidis et al. 2003). The interaction with this shell was identified as the main reason why metal-carrying rising Rayleigh–Taylor structures were unable to penetrate into the H-shell in the 2D simulations of Kifonidis et al. (2003). Here we will integrate the equation of motion for spherical clumps decelerated by the drag force in the 3D environment (in contrast to the 3D case of a spherical clump, the equation of motion of a 2D toroidal object defies an analytic solution due to the explicit dependence of the clump volume and cross-sectional area on the radial coordinate position r).

We assume that the He layer is homogeneous but has a time-dependent density ρ(t). This is an acceptable simplification because numerical simulations show that the radial density variations in this layer are much smaller than the time-dependent decrease of the density due to the expansion of the star during the first three hours of the explosion. In the 3D case this leaves us just with a time dependence of the drag force, F_drag = F_drag(t), and thus

$$\begin{equation} m {{{d}}v(r)\over {{d}}t}= -{1\over 2} \, C_\mathrm{s} A_\mathrm{s}(t)\rho (t) v^2(t) \,. \end{equation} \tag{ 8 }$$

Again we consider m = const and C_s = const and v(t) as the clump velocity relative to the surrounding medium. We further assume that the density of the helium layer evolves like ρ(t) = ρ₀(t₀/t)^k and the pressure and density are linked by P(t) = Kρ^γ(t). Using the same pressure–density relation for the medium inside the clumps, P_s = K_sρ^γ_s, but with a different constant K_s ≠ K (corresponding to an entropy of the clumps different from that of the surroundings), the request of pressure equilibrium between clumps and environment allows us to calculate the time evolution of the clump radius,

$$\begin{equation} r_\mathrm{s}(t) = r_{\mathrm{s},0} \left({K_\mathrm{s}\over K}\right)^{\! {1\over 3\gamma }} \left({\rho _{\mathrm{s},0}\over \rho _0 }\right)^{\! {1\over 3}} \left({t\over t_0}\right)^{\! {k\over 3}}, \end{equation} \tag{ 9 }$$

and therefore the time-dependent cross-sectional area that the clumps present to the surrounding medium: A_s(t) = A_s,0(r/r_s,0)² ≡ A*(t/t₀)^2k/3 with constant A*. Inserting this and ρ(t) into Equation (8), we can integrate the resulting differential equation by separation of variables. This yields the time evolution of the velocity (for k ≠ 3) as

$$\begin{equation} v(t)= v_0 \,\left\lbrace 1\,+\, {C_\mathrm{s}A^\ast \over 2m} v_0t_0\rho _0 {1\over 1-k/3} \left[ \left({t\over t_0}\right)^{\! 1-k/3} - 1 \right] \right\rbrace ^{-1}. \end{equation} \tag{ 10 }$$

We specialize now to k = 2 and $\gamma = {4\over 3}$, which is in reasonable agreement with the simulations (Kifonidis et al. 2003). Rewriting the constant factor of the second term in the denominator by using the expressions for clump area and volume, we end up with

$$\begin{equation} v(t) = v_0 \,\left\lbrace 1\,+\, {9\over 8}\left(\! {4\pi \over 3}\! \right)^{\!\! {1\over 3}}\!\! \left({K_{\mathrm s}\over K}\right)^{\!\! {1\over 2}} \! C_\mathrm{s} v_0t_0\left({\rho _0\over m}\right)^{\!\! {1\over 3}} \!\! \left[ \left({t\over t_0}\right)^{\!\! {1\over 3}}\!\! - 1 \right] \right\rbrace ^{-1} \!\!\!. \end{equation} \tag{ 11 }$$

The fastest clumps, which are also the most massive ones with m ∼ 4 × 10⁻³ M_☉ in our 3D simulations, enter the dense helium shell with a velocity of v₀ ∼ 2000–4000 km s⁻¹ (relative to the ambient medium) at t₀ ∼ 300 s. At this time the density of the medium in the shell is ρ₀ ∼ 0.1 g cm⁻³. The clumps typically have roughly three times lower entropy than the medium in the helium shell so that K_s/K ∼ 1/3 (most of the numerical values used here can be verified by closely inspecting the plots in Kifonidis et al. 2003). Adopting for the drag coefficient plausible numbers⁵ of C_s ∼ 0.05 ... 0.2, we obtain at t = 9000 s velocities of v ∼ 1200–3100 km s⁻¹ for the clump motion through the helium material, corresponding to a deceleration by roughly a factor of 2. This means that the fastest objects are expected to still speed with ∼2400–4400 km s⁻¹ relative to an observer in the lab frame. This is in good agreement with our hydrodynamical simulations and can thus be interpreted as support of our picture of the clump propagation under the influence of the drag force.⁶ In contrast, in 2D models the maximum clump velocities are reduced to slightly more than 1500 km s⁻¹ in the lab frame (see Figures 6, 8 and Kifonidis et al. 2003), which is equal to the expansion velocity of the He-shell matter so that the clumps get stuck in their surroundings and are unable to penetrate into the hydrogen envelope (Figure 7).