IS STRONG SASI ACTIVITY THE KEY TO SUCCESSFUL NEUTRINO-DRIVEN SUPERNOVA EXPLOSIONS?

Recent simulations in two dimensions with sophisticated neutrino transport have demonstrated that the neutrino-driven mechanism, supported by hydrodynamic instabilities in the postshock layer, can yield supernova explosions at least for some progenitor stars (11.2 and 15 M_☉ ones in Marek & Janka 2009; Müller et al. 2012). The explosions occur relatively late after bounce and tend to have fairly low energy, being "marginal" or only slightly above the "critical threshold" in this sense. Suwa et al. (2010) obtained a similar explosion for a 13 M_☉ progenitor in their axisymmetric simulations. However, the Oak Ridge group has found stronger and earlier explosions for a wider range of progenitors (Bruenn et al. 2009), while in purely Newtonian simulations with multi-dimensional neutrino diffusion (including energy dependence but without energy-bin coupling) Burrows et al. (2006, 2007) could not see any success of the delayed neutrino-driven mechanism.

While the reason for the discrepant results of these simulations cannot be satisfactorily understood on the basis of published results, the marginality of the two-dimensional (2D) explosions of the Garching group and the lack of neutrino-driven explosions in the simulations by Burrows et al. (2006, 2007) raises an important question about the influence of the third spatial dimension on the post-bounce evolution of collapsing stellar cores. Three-dimensional (3D) fluid dynamics with their inverse turbulent energy cascade compared to the 2D case are likely to change the flow pattern on large scales as well as small scales. They could have an influence on the existence and the growth rate of non-radial hydrodynamic instabilities in the supernova core even in the absence of stellar rotation (see, e.g., Iwakami et al. 2008), but in particular with a moderate amount of angular momentum in the progenitor star (e.g., Blondin & Mezzacappa 2006a; Iwakami et al. 2009; Fernández 2010), and thus could lead to differences in the hydrodynamic and thermodynamic conditions for the operation of the neutrino-heating mechanism. In particular, 3D flows might cause important changes of the dwell time of postshock matter in the layer where neutrinos deposit energy, which is a crucial aspect deciding about the viability and efficiency of the neutrino-driven supernova mechanism (some aspects of this were discussed by Murphy & Burrows 2008; Marek & Janka 2009).

Indeed, employing a simplified treatment of neutrino effects by including local neutrino-cooling and heating terms for a chosen value of the neutrino luminosity and spectral temperature instead of solving the computationally intense neutrino transport, Nordhaus et al. (2010) found considerably easier and earlier explosions in three dimensions than in two dimensions. In the context of the concept of a critical value of the neutrino luminosity that (for a given mass accretion rate onto the stalled supernova shock) has to be exceeded to obtain neutrino-driven explosions (Burrows & Goshy 1993; Janka & Müller 1996; Janka 2001; Yamasaki & Yamada 2005; Murphy & Burrows 2008; Pejcha & Thompson 2011; Fernández 2012), they quantified the improvement of three dimensions relative to two dimensions by a 15%–25% reduction of the critical luminosity value. In particular, they observed that 3D postshock convection leads to higher average entropies in the neutrino-heating layer, thus improving the conditions for shock revival due to a significant stretching of the residence time of matter in the layer where it gains energy from neutrinos. Very recently, Takiwaki et al. (2011) reported enhanced maximum dwell times of a small fraction of the material in the gain region in a 3D simulation compared to the 2D case of an 11.2 M_☉ star, but could not unambiguously link this effect to an easier explosion of the 3D model. In particular, their 3D simulation showed a shock expansion that was more delayed than in the 2D run, and the 3D conditions did not appear more favorable for an explosion with respect to a variety of quantities like the net heating rate, the heating timescale, or the profiles of maximum and minimum entropies.

In this paper, we present a comparative investigation for 11.2 M_☉ and 15 M_☉ progenitors in one, two, and three dimensions along the lines of the study by Nordhaus et al. (2010), varying the driving neutrino luminosities used in time-dependent collapse simulations of the two stars. While our results for spherically symmetric (1D) and 2D models basically confirm the dimension-dependent differences found by Murphy & Burrows (2008) and Nordhaus et al. (2010), our calculations neither exhibit a strict 1D–2D–3D hierarchy of the average entropy in the gain layer, nor do they show any clear signs that 3D effects facilitate the development of the explosion better than non-radial motions in 2D. Attempting to understand the reason for this puzzling finding we vary the resolution of the spherical coordinate grid used for our 2D and 3D simulations. The outcome of these studies reflects the action of the energy flow within the turbulent energy cascade. The latter transports the driving energy provided by neutrino heating and gravitational energy release in the accretion flow from small to large scales in two dimensions and opposite in three dimensions. Models in two dimensions show growing large-scale asymmetry and quasi-periodic time variability and clearly explode easier with higher resolution, whereas in three dimensions better resolved models are observed to become more similar to the 1D case and thus to be farther away from an explosion. This suggests that the success of the neutrino-driven mechanism could be tightly linked to the initiation of strong non-radial mass motions in the neutrino-heated postshock layer on the largest possible scales, implying that the easier explosions of our 2D models with higher resolution are a consequence of more violent activity due to the standing accretion-shock instability (SASI; Blondin et al. 2003), whereas the better resolved 3D models for the employed artificial setup of supernova-core conditions tend to reveal considerably reduced amplitudes of low-order spherical-harmonics modes of non-radial deformation and thus behave more similar to the 1D case.

The paper is structured as follows. In Section 2, we briefly describe our numerics and implementation of neutrino source terms. In Section 3, we give an overview of the simulations presented in this paper. Our investigations of the dependence of the critical luminosity on the spatial dimension are presented in Section 4 and results of resolution studies are presented in Section 5. An interpretation of our findings will follow in Section 6. Section 7 contains the summary and conclusions. In the Appendix, we present 1D simulations that document our efforts to reproduce the results of previous, similar studies in the literature by straightforwardly applying the neutrino treatment described in these works.

We solve the equations of hydrodynamics reflecting the conservation of mass, momentum, and energy,

$$\begin{equation} \frac{\partial \rho }{\partial t}+\mathbf {\nabla }\cdot (\rho \mathbf {v})=0\,, \end{equation} \tag{ 1 }$$

$$\begin{equation} \frac{\partial \rho \mathbf {v}}{\partial t}+ \mathbf {\nabla } \cdot (\rho \mathbf {v} \otimes \mathbf {v}) + \mathbf {\nabla } P = -\rho \mathbf {\nabla } \Phi \,, \end{equation} \tag{ 2 }$$

$$\begin{equation} \frac{\partial e}{\partial t}+ \mathbf {\nabla } \cdot \left[\left(e + P\right) \mathbf {v} \right] = -\rho \mathbf {v} \cdot \mathbf {\nabla } \Phi + \rho \left(Q_\nu ^+ - Q_\nu ^-\right) \,, \end{equation} \tag{ 3 }$$

where ρ is the mass density, $\mathbf {v}$ is the fluid velocity, Φ is the gravitational potential, P is the pressure, and e is the total (internal+kinetic) fluid energy density. These equations are integrated in a conservative form (for which reason the energy equation is solved for the total energy density) using the explicit, finite-volume, higher-order Eulerian, multi-fluid hydrodynamics code Prometheus (Fryxell et al. 1991; Müller et al. 1991a, 1991b). It is a direct implementation of the piecewise parabolic method (PPM) of Colella & Woodward (1984) using the Riemann solver for real gases developed by Colella & Glaz (1985) and the directional splitting approach of Strang (1968) to treat the multi-dimensional problem. In order to prevent odd–even coupling (Quirk 1994) we switch from the original PPM method to the HLLE solver of Einfeldt (1988) in the vicinity of strong shocks. The advection of nuclear species is treated by the consistent multi-fluid advection scheme of Plewa & Müller (1999).

To facilitate comparison with Nordhaus et al. (2010) we also employ the high-density equation of state (EoS) of Shen et al. (1998) and do not include general relativistic corrections. We use the monopole approximation of the Poisson equation to treat Newtonian self-gravity.

To make our extensive parameter study possible, we use the local source terms applied by Murphy & Burrows (2008) and Nordhaus et al. (2010) instead of detailed neutrino transport (see Janka 2001 for a derivation of these source terms). In this approach, the neutrino-heating and cooling rates Q⁺_ν and Q⁻_ν are given by

$$\begin{eqnarray} Q^{+}_{\nu } & = & 1.544\times 10^{20} \left(\frac{L_{\nu _{e}}}{10^{52} \,{\rm erg s}^{-1}}\right) \left(\frac{T_{\nu _{e}}}{4\,{\rm MeV}}\right)^{2} \\ & & \times \left(\frac{100\,{\rm km}}{r}\right)^2 (Y_{n}+Y_{p}) \, e^{-\tau _{{\rm eff}}} \left[\frac{{\rm erg}}{{\rm g s}}\right] \,, \nonumber \end{eqnarray} \tag{ 4 }$$

$$\begin{equation} Q^{-}_{\nu } = 1.399\times 10^{20} \left(\frac{T}{2\,{\rm MeV}}\right)^{6} (Y_{n}+Y_{p}) \, e^{-\tau _{{\rm eff}}} \left[\frac{{\rm erg}}{{\rm g s}}\right] \,. \end{equation} \tag{ 5 }$$

These approximations depend on local quantities, namely, the density ρ, the temperature T, the distance from the center of the star r, and the neutron and proton number fractions Y_n and Y_p, respectively. In Equation (4), the electron–neutrino luminosity $L_{\nu _{e}}$ is a parameter and is assumed to be equal to the electron–antineutrino luminosity $L_{\bar{\nu }_{e}} = L_{\nu _{e}}$. The neutrino temperature $T_{\nu _{e}}$ is set to 4 MeV.

The employed source terms, Equations (4) and (5), without the factors $e^{-\tau _{{\rm eff}}}$ are valid for optically thin regions only. In order to model the transition to neutrino trapping at high optical depths, we follow Nordhaus et al. (2010) and multiply the heating and cooling terms by $e^{-\tau _{{\rm eff}}}$ to suppress them in the inner opaque core of the proto-neutron star. Here, the optical depth for electron neutrinos and antineutrinos is defined as

$$\begin{equation} \tau _{{\rm eff}}(r) = \int _{r}^{\infty } \kappa _{{\rm eff}}(r^{\prime }) \, {d}r^{\prime } . \end{equation} \tag{ 6 }$$

The effective opacity κ_eff was derived by Janka (2001) and given in Murphy et al. (2009):

$$\begin{eqnarray} \kappa _{{\rm eff}} &\approx& 1.5\times 10^{-7}\cdot X_{n,p} \left(\frac{\rho }{10^{10}\,{\rm g cm}^{-3}}\right)\left(\frac{T_{\nu _{e}}}{4 \,{\rm MeV}}\right)^{2} {\rm cm}^{-1},\nonumber\\ \end{eqnarray} \tag{ 7 }$$

where X_{n, p} = (1/2)(Y_n + Y_p) accounts for composition averaging. In multi-dimensional simulations, we evaluate the radial integration for the optical depth τ_eff independently for each latitude θ in two dimensions and each pair of latitudinal and azimuthal angles (θ,ϕ) in three dimensions. Note that in Murphy & Burrows (2008) the exponential suppression factor is absent in the heating and cooling terms (or was not mentioned), which otherwise agree with ours, while no definition of the effective optical depth τ_eff is given in Nordhaus et al. (2010). The factor (Y_n + Y_p) in Equations (4) and (5) is included in Nordhaus et al. (2010), but not in Murphy & Burrows (2008) or Murphy et al. (2009).

The time period from the onset of the collapse until 15 ms after bounce was tracked with the Prometheus-Vertex code (Rampp & Janka 2002) in one dimension including its detailed, multi-group neutrino transport, all relevant neutrino reactions, and a "flashing treatment" for an approximative description of nuclear burning during infall. This means that until 15 ms after bounce we describe neutrino effects including the evolution of the electron fraction Y_e with high sophistication. At 15 ms after core bounce we switch to the simple neutrino-heating and cooling terms and upon mapping from one dimension to two dimensions also impose (on the whole computational grid) random zone-to-zone seed perturbations with an amplitude of 1% of the density to break spherical symmetry.

Although during the subsequent evolution we apply the heating and cooling expressions of Equations (4) and (5) following Nordhaus et al. (2010) and Murphy & Burrows (2008), we refrain from adopting their treatment of changes of the electron fraction Y_e. Following a suggestion by Liebendörfer (2005), they prescribed Y_e simply as a function of density, Y_e(ρ), instead of solving a rate equation with source terms for electron neutrino and antineutrino production and destruction consistently with the expressions employed for neutrino heating and cooling. Liebendörfer (2005) found that such a parameterization, supplemented by a corresponding entropy source term (and a neutrino pressure term in the equation of motion), yields results in good agreement with 1D simulations including neutrino transport during the collapse phase until the moment of bounce-shock formation. A universal Y_e–ρ relation, which serves as the basis of this approximation and can be inferred for the infalling matter during the homologous collapse phase, however, applies neither to the evolution of the shocked accretion flow in the post-bounce phase nor to expanding neutrino-heated gas (see, e.g., Figure 4.9 in Thielemann et al. 2011). For example, a comparison with supernova models with detailed neutrino transport shows that the Y_e–ρ relation fitted to the homologous phase overestimates the deleptonization of the gas in most of the gain layer but underestimates the lepton loss of matter in the cooling layer and neutrinospheric region. Moreover, the question arises: How the lepton evolution shall be treated in matter that reexpands and thus moves from high to low density? Even more, the entropy source term introduced in Equation (5) of Liebendörfer (2005) is designed to specifically account for gas entropy changes due to neutrino production by electron captures and subsequent energy transfers in (neutrino-electron) scatterings. The corresponding energy loss or gain rate of the medium through the escape or capture of electron neutrinos with mean energy $E_{\nu _e}$, which is given by $\delta Q_{\nu _e}/\delta t = E_{\nu _e}\,\delta Y_e/\delta t$, is not included in the heating and cooling terms Q⁺_ν and Q⁻_ν of Equations (4) and (5) of the present work. Adding it as an extra term would imply partial double counting of the energy exchange via electron neutrinos, and omitting it means overestimating the entropy increase in infalling, deleptonizing matter and underestimating entropy gains of decompressed gas with growing Y_e. Because of the long list of such inconsistencies, whose implications are hard to judge or control, we decided to ignore Y_e changes of the stellar medium in our post-bounce simulations (except for the models discussed in the Appendix).

This choice can be justified, but it is certainly not a perfect approach because it may also exclude important effects in real supernova cores, whose physical processes require the inclusion of neutrino transport for an accurate description of the energy and lepton-number evolution. One of the undesired consequences of keeping Y_e fixed in the accretion flow after core bounce is an overestimation of the electron pressure in the gas settling onto the forming neutron star. In order to enforce more compression and thus to ensure close similarity of our results to the 1D models studied by Nordhaus et al. (2010) and Murphy & Burrows (2008), we have to enhance the net cooling of the accreted matter by reducing the effective opacity κ_eff by a factor of 2.7 compared to the value given in Murphy et al. (2009) and in Equation (7). This reduction factor is chosen such that our simulations reproduce the minimum values of the critical luminosity found to be necessary for triggering explosions of the 15 M_☉ progenitor (2.6 × 10⁵² erg s⁻¹) and for the 11.2 M_☉ star (1.3 × 10⁵² erg s⁻¹) in the 1D simulations of Murphy & Burrows (2008). Without the reduction factor of κ_eff, our models turn out to explode too easily because of weak cooling (see the results in the Appendix). We stress that any exponential suppression factor of the heating and cooling rates in Equations (4) and (5) is a pragmatic and ad hoc procedure to bridge the transition from the optically thin to the optically thick regime, where neutrino transport is most complicated. From transport theory neither the exponential factor nor the exact definition of the optical depth of the exponent can be rigorously derived.

In our reference set of standard simulations, we employ 400 non-equidistant radial zones, which are distributed from the center to an outer boundary at 9000 km. The latter is sufficiently far out to ensure that the gas there remains at rest during the simulated evolution periods. The radial zones are chosen such that the resolution Δr/r is typically better than 2%. In the multi-dimensional models with standard resolution, we employ a polar grid with an angular bin size of 3° (60 θ- and 120 ϕ-zones). For the high-resolution 1D and multi-dimensional models of Sections 4 and 5 we compute with up to 800 radial zones and in two dimensions with an angular resolution down to 05, and in three dimensions down to 15. The additional radial grid zones are distributed such that the region between 20 and 400 km, i.e., not only the cooling layer around the neutrinosphere but also the gain layer between gain radius and shock, are resolved significantly better. While full convergence of the 1D results requires 600 radial zones or more, most of the plots comparing models for different dimensions show simulations with our standard resolution of only 400 radial cells (unless stated otherwise) because this is the limit for which we could perform a larger set of 3D runs in acceptable time. To avoid an extremely restrictive Courant–Friedrich–Lewy time step in our 3D calculations, we simulate the inner core above a density of ρ = 10¹² g cm⁻³ in spherical symmetry.

Our models are based on the 15 M_☉ progenitor star s15s7b2 of Woosley & Weaver (1995) and an 11.2 M_☉ progenitor of Woosley et al. (2002). The calculations for these progenitors with our standard angular resolution of 3° are summarized in Table 1. This table is arranged such that horizontal rows have the same driving luminosities for simulations performed in different dimensions. Varying the prescribed driving luminosity $L_{\nu _e}$ from run to run we present for each of the 11.2 M_☉ and 15 M_☉ progenitors several 1D, 2D, and 3D simulations. All of the 1D and 2D simulations cover at least 1 s after bounce. The nonexploding 3D simulations with standard angular resolution of 3° were not stopped until at least 600 ms after bounce. For models with higher resolution the simulation times are given in Table 2.

Table 1. 11.2 M_☉ and 15 M_☉ Results with Standard Grid of 400 Radial Zones and 3° Angular Resolution

	1D		2D		3D
$L_{\nu _e}$a	texpb	$\dot{M}_{{\rm exp}}$c	texp	$\dot{M}_{{\rm exp}}$	texp	$\dot{M}_{{\rm exp}}$
(1052 erg s−1)	(ms)	(M☉ s−1)	(ms)	(M☉ s−1)	(ms)	(M☉ s−1)
11.2 M☉
0.8	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅
0.9	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	731	0.085
1.0	⋅⋅⋅	⋅⋅⋅	563	0.082	537	0.086
1.1	⋅⋅⋅	⋅⋅⋅	461	0.091
1.2	⋅⋅⋅	⋅⋅⋅	357	0.104	319	0.112
1.3	819	0.082	307	0.114
1.4	499	0.088	241	0.126	232	0.130
1.5	380	0.100	232	0.130
1.6	345	0.106	203	0.137
15 M☉
2.0	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅
2.1	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅
2.2	⋅⋅⋅	⋅⋅⋅	876	0.197	612	0.226
2.3	⋅⋅⋅	⋅⋅⋅	428	0.261	426	0.261
2.4	⋅⋅⋅	⋅⋅⋅	442	0.255
2.5	⋅⋅⋅	⋅⋅⋅	283	0.313	281	0.314
2.6	710	0.215	285	0.312
2.7	489	0.247	262	0.316
2.8	390	0.271	242	0.322	236	0.324
2.9	281	0.314	235	0.325
3.0	258	0.318	236	0.324
3.1	248	0.320	220	0.327

Notes. ^aElectron–neutrino luminosity. ^bTime after bounce of onset of explosion. A "−" symbol indicates that the model does not explode during the simulated period of evolution. ^cMass accretion rate at onset of explosion.

Download table as:  ASCIITypeset image

Table 2. Multi-dimensional Models with Different Resolution

Massa	Dimb	$L_{\nu _e}$c	Ang.d	Nre	texpf	$\dot{M}_{{\rm exp}}$g	tsimh
(M☉)		(1052 erg s−1)	Res.		(ms)		(ms)
11.2	2D	0.8	3°	400	⋅⋅⋅	⋅⋅⋅	1017
11.2	2D	0.8	2°	400	⋅⋅⋅	⋅⋅⋅	979
11.2	3D	0.8	3°	400	⋅⋅⋅	⋅⋅⋅	915
11.2	3D	0.8	2°	400	⋅⋅⋅	⋅⋅⋅	758
11.2	2D	0.9	3°	400	⋅⋅⋅	⋅⋅⋅	1006
11.2	2D	0.9	2°	400	⋅⋅⋅	⋅⋅⋅	985
11.2	3D	0.9	3°	400	731	0.085	954
11.2	3D	0.9	2°	400	⋅⋅⋅	⋅⋅⋅	819
11.2	2D	1.0	3°	400	563	0.082	1053
11.2	2D	1.0	2°	400	527	0.086	1053
11.2	3D	1.0	3°	400	537	0.086	684
11.2	3D	1.0	2°	400	572	0.082	761
15	2D	2.0	3°	400	⋅⋅⋅	⋅⋅⋅	1016
15	2D	2.0	2°	400	⋅⋅⋅	⋅⋅⋅	829
15	2D	2.0	15	400	⋅⋅⋅	⋅⋅⋅	1016
15	2D	2.0	1°	400	⋅⋅⋅	⋅⋅⋅	1016
15	3D	2.0	3°	400	⋅⋅⋅	⋅⋅⋅	723
15	3D	2.0	2°	400	⋅⋅⋅	⋅⋅⋅	524
15	2D	2.1	3°	400	⋅⋅⋅	⋅⋅⋅	1016
15	2D	2.1	2°	400	⋅⋅⋅	⋅⋅⋅	829
15	2D	2.1	15	400	719	0.210	1016
15	2D	2.1	1°	400	575	0.232	1016
15	2D	2.1	05	400	657	0.220	1016
15	3D	2.1	3°	400	⋅⋅⋅	⋅⋅⋅	767
15	3D	2.1	2°	400	⋅⋅⋅	⋅⋅⋅	815
15	3D	2.1	15	400	⋅⋅⋅	⋅⋅⋅	777
15	2D	2.2	3°	400	876	0.197	1016
15	2D	2.2	2°	400	557	0.238	825
15	2D	2.2	15	400	556	0.239	1016
15	2D	2.2	1°	400	424	0.262	1016
15	2D	2.2	05	400	365	0.288	1016
15	3D	2.2	3°	400	612	0.226	932
15	3D	2.2	2°	400	⋅⋅⋅	⋅⋅⋅	925
15	2D	2.1	3°	600	⋅⋅⋅	⋅⋅⋅	1016
15	2D	2.1	3°	800	⋅⋅⋅	⋅⋅⋅	1016
15	2D	2.1	2°	600	⋅⋅⋅	⋅⋅⋅	1016
15	2D	2.1	2°	800	⋅⋅⋅	⋅⋅⋅	1016
15	2D	2.1	15	600	⋅⋅⋅	⋅⋅⋅	1016
15	2D	2.1	15	800	⋅⋅⋅	⋅⋅⋅	1016
15	2D	2.1	1°	600	780	0.203	1016
15	2D	2.1	1°	800	⋅⋅⋅	⋅⋅⋅	1016
15	2D	2.1	05	600	749	0.211	1016
15	2D	2.1	05	800	886	0.198	1016
15	3D	2.1	3°	600	⋅⋅⋅	⋅⋅⋅	946
15	3D	2.1	3°	800	⋅⋅⋅	⋅⋅⋅	958
15	2D	2.2	3°	600	⋅⋅⋅	⋅⋅⋅	1016
15	2D	2.2	3°	800	⋅⋅⋅	⋅⋅⋅	1016
15	2D	2.2	2°	600	683	0.218	1016
15	2D	2.2	2°	800	⋅⋅⋅	⋅⋅⋅	1016
15	2D	2.2	15	600	713	0.215	1016
15	2D	2.2	15	800	863	0.196	1016
15	2D	2.2	1°	600	761	0.214	1016
15	2D	2.2	1°	800	961	0.194	1016
15	2D	2.2	05	600	588	0.232	1016
15	2D	2.2	05	800	643	0.224	1016
15	3D	2.2	3°	600	803	0.202	975
15	3D	2.2	3°	800	857	0.195	977
15	3D	2.2	2°	600	⋅⋅⋅	⋅⋅⋅	988
15	3D	2.2	15	800	⋅⋅⋅	⋅⋅⋅	898

Notes. 3D models are indicated in bold. ^aProgenitor model. ^bDimensionality. ^cElectron–neutrino luminosity. ^dAngular resolution. ^eNumber of radial zones. ^fTime of onset of explosion. ^gMass accretion rate at onset of explosion. ^hSimulation time.

Download table as:  ASCIITypeset image

In Table 1 the time of the onset of an explosion, t_exp, and the mass accretion rate at that time, $\dot{M}_\mathrm{exp}$, are listed as characteristic quantities of the models. The beginning of the explosion is defined as the moment of time t_exp when the shock reaches an average radius of 400 km (and does not return later on), while nonexploding models are denoted by an "−" symbol. In multi-dimensional simulations, the corresponding shock position is defined as the surface average over all angular directions, 〈R_S〉 ≡ (1/4π)∮dΩ R_S(Ω). The lowest driving luminosity yielding an explosion for a given value of the mass accretion rate is termed the critical luminosity (Burrows & Goshy 1993; Murphy & Burrows 2008). We determine the mass accretion rate $\dot{M}(r) = 4\pi r^2\rho (r) \left|v(r)\right|$ at the time of the onset of the explosion just exterior to the shock, i.e., at a radius of 500 km, where the infalling envelope is spherical (except for the small seed perturbations imposed on the multi-dimensional models). In Figure 1 the mass accretion rates of the 11.2 M_☉ and 15 M_☉ progenitors are depicted for the nonexploding runs with the lowest driving neutrino luminosities. Because the shock can be largely deformed in multi-dimensional simulations and its outermost parts can extend beyond 500 km (and thus impede a clean determination of the mass accretion rate) when its average radius just begins to exceed 400 km, we refer to the functions plotted in Figure 1 for defining the mass accretion rate at the time when the explosion sets in.

Zoom In Zoom Out Reset image size

Based on steady-state solutions of neutrino-heated and cooled accretion flows between the stalled shock and the proto-neutron star surface, Burrows & Goshy (1993) identified a critical condition that can be considered to separate exploding from nonexploding models. They found that for a given value of the mass infall rate $\dot{M}$ onto the accretion shock steady-state solutions cannot exist when the neutrino-heating rate in the gain layer is sufficiently large, i.e., for neutrino luminosities above a threshold value L_ν. This result can be coined in terms of a critical condition $L_\nu (\dot{M})$ expressing the fact that either the driving luminosity has to be sufficiently high or the damping mass accretion rate enough low for an explosion to become possible. The critical $\dot{M}$–L_ν curve was interpreted by Burrows & Goshy (1993) as a separating line between the region above, where due to the non-existence of steady-state accretion solutions explosions are expected to take place, from the region below, where neutrino-energy input behind the shock is insufficient to accelerate the stalled shock outward and thus to trigger an explosion.

This interpretation of the steady-state results was consistent with hydrodynamical simulations of collapsing and exploding stars in one dimension and two dimensions by Janka & Müller (1996). Performing a more extended parameter study than the latter work, Murphy & Burrows (2008) explored the concept of a critical condition systematically with time-dependent hydrodynamical models. They showed that a critical luminosity indeed separates explosion from accretion and confirmed that this value is lowered by ∼30% when going from spherical symmetry to two dimensions, at least when a fixed driving luminosity is adopted and feedback effects of the hydrodynamics on the neutrino emission are ignored. Some 2D effects including possible consequences of rotation were discussed before on grounds of steady-state models by Yamasaki & Yamada (2005, 2006), while Janka (2001) tried to include time-dependent aspects of the shock-revival problem and took into account an accretion component of the neutrino luminosity in addition to the fixed core component. The influence of such an accretion contribution was more recently also estimated by Pejcha & Thompson (2011), who solved the 1D steady-state accretion problem between the neutron star and the accretion shock along the lines of Burrows & Goshy (1993), but attempted to obtain a deeper understanding of the critical condition by comprehensively analyzing the structure of the accretion layer and of the limiting steady-state solution in dependence of the stellar conditions. They found that the critical value for the neutrino luminosity is linked to an "antesonic condition" in which the ratio of the adiabatic sound speed to the local escape velocity in the postshock layer reaches a critical value above which steady-state solutions of neutrino-heated accretion flows cannot be obtained. By performing high-resolution hydrodynamic simulations, Fernández (2012) found that radial instability is a sufficient condition for runaway expansion of an initially stalled core-collapse supernova shock if the neutrinospheric parameters do not vary with time and if heating by the accretion luminosity is neglected. However, the threshold neutrino luminosities for the transition to runaway instability are in general different from the limiting values for steady-state solutions of the kind discussed by Burrows & Goshy (1993) and Pejcha & Thompson (2011). Nordhaus et al. (2010) generalized the hydrodynamic investigations of Murphy & Burrows (2008) to include 3D models and found another reduction of the threshold luminosity for explosion by 15%–25% compared to the 2D case.

Despite the basic agreement of the outcome of these investigations it should be kept in mind that it is not ultimately clear whether the simple concept of a critical threshold condition separating explosions from failures (and the dependences of this threshold on dimension and rotation, for example) holds beyond the highly idealized setups considered in the mentioned works. None of the mentioned systematic studies by steady-state or hydrodynamic models was able to include adequately the complexity of the feedback between hydrodynamics and neutrino transport physics. In particular, none of these studies could yield the proof that the non-existence of a steady-state accretion solution for a given combination of mass accretion rate and neutrino luminosity is equivalent to the onset of an explosion. The latter requires the persistence of sufficiently strong energy input by neutrino heating for a sufficiently long period of time. This is especially important because Pejcha & Thompson (2011) showed that the total energy in the gain layer is still negative even in the case of the limiting accretion solution that corresponds to the critical luminosity.¹ Within the framework of simplified modeling setups, however, the question cannot be answered whether such a persistent energy input can be maintained in the environment of the supernova core.

Following the previous investigations by Murphy & Burrows (2008) and Nordhaus et al. (2010) we performed hydrodynamical simulations that track the post-bounce evolution of collapsing stars for different, fixed values of the driving neutrino luminosity. Since the mass accretion rate decreases with time according to the density profile that is characteristic of the initial structure of the progenitor core (see Figure 1 for the 11.2 and 15 M_☉ stars considered in this work), each model run probes the critical value of $\dot{M}_{\mathrm{exp}}$ at which the explosion becomes possible for the chosen value of $L_\nu = L_{\nu _e} = L_{\bar{\nu }_e}$. The collection of value pairs ($\dot{M}_{\mathrm{exp}}$, $L_{\nu _e}$) defines a critical curve $L_{\nu }(\dot{M})$. These are shown for our 1D, 2D, and 3D studies with standard resolution for both progenitor stars in the left panel of Figure 2 and in the case of the 15 M_☉ star can be directly compared with Figure 1 of Nordhaus et al. (2010). Table 1 lists, as a function of the chosen $L_{\nu _e}$, the corresponding times t_exp when the onset of the explosion takes place and the mass accretion rate has the value of $\dot{M}_{\mathrm{exp}}$. The post-bounce evolution of a collapsing star proceeds from high to low mass accretion rate (Figure 1), i.e., from right to left on the horizontal axis of the left panel of Figure 2. When $\dot{M}$ reaches the critical value for the given $L_{\nu _e}$, the model develops an explosion. The right panel of Figure 2 visualizes the functional relations between the neutrino luminosities $L_{\nu _e}$ and the explosion times t_exp for both progenitors and for the simulations with different dimensions.

Zoom In Zoom Out Reset image size

At first glance Figure 2 reproduces basic trends that are visible in Figure 17 of Murphy & Burrows (2008) and in Figure 1 of Nordhaus et al. (2010). For example, the critical luminosity increases for higher mass accretion rate and the values for spherically symmetric models are clearly higher than those for 2D simulations. However, a closer inspection reveals interesting differences compared to the previous works.

• 1.  
  In general the slopes of our critical $L_{\nu _e}(\dot{M})$ relations appear to be considerably steeper and in the case of the 15 M_☉ star they exhibit a very steep rise above $\dot{M}\approx 0.31\,M_\odot$ s⁻¹. This means that explosions for higher mass accretion rates are much harder to obtain and therefore the tested driving luminosities in our simulations do not lead to explosions earlier than about 200 ms after bounce, independent of whether the modeling is performed in one dimension, two dimensions, or three dimensions. In contrast, the critical curves given by Nordhaus et al. (2010) show a moderately steep increase over the whole range of plotted mass accretion rates between about 0.1 M_☉ s⁻¹ and more than 0.5 M_☉ s⁻¹.
• 2.  
  Non-radial flows in the 2D case, by which the residence time of accreted matter in the gain layer could be extended or more matter could be kept in the gain region, reduce the critical luminosities by at most ∼15% of the 1D values for the 15 M_☉ star and ≲25% for the 11.2 M_☉ model, which is a somewhat smaller difference than found by Murphy & Burrows (2008) and Nordhaus et al. (2010).
• 3.  
  Most important, however, is the fact that we cannot confirm the observation by Nordhaus et al. (2010) that the 3D case provides considerably more favorable conditions for an explosion than the 2D case. Our critical curves for the 2D and 3D cases nearly lie on top of each other. There are minor improvements of the explosion conditions in three dimensions visible in both panels of Figure 2 and the numbers of Table 1, e.g., a 10% reduction of the smallest value of L_ν for which an explosion can be obtained for the 11.2 M_☉ star, a ∼260 ms earlier explosion for the lowest luminosity driving the 15 M_☉ explosion (2.2 × 10⁵² erg s⁻¹), and a tendency of slightly faster 3D explosions for all tested luminosities (see Figure 3 and Table 1). All of these more optimistic 3D features, however, will disappear for simulations with higher resolution as we will see in Section 5.

Before we discuss the origin of the differences between our results and those of Murphy & Burrows (2008) and Nordhaus et al. (2010) we would like to remark that the kinks and even nonmonotonic parts of the curves shown in Figure 2 in particular for the multi-dimensional cases are connected to our definition of the explosion time as being the moment when the mean shock radius exceeds 400 km. Especially in cases where the shock deformation is large (which is an issue mainly in some of the 2D simulations) this definition is associated with significant uncertainty in the determination of the exact explosion time t_exp and therefore also of the corresponding mass accretion rate $\dot{M}_{\mathrm{exp}}$.

In addition to the results for our standard resolution, Figure 2 presents 1D models with higher radial resolution (600 and 800 radial zones). The critical curves with better zoning exhibit the same overall slopes as those of the standard runs, but there is a slight shift toward higher values of the critical luminosity (or, equivalently, a small shift to lower $\dot{M}_{\mathrm{exp}}$ and later t_exp). This is caused by somewhat larger neutrino-energy losses from the cooling layer with better radial resolution, an effect which makes explosions more difficult. We stress that this resolution-dependent cooling effect is a consequence of the simplified neutrino-loss treatment. The employed cooling rate is not able to reproduce real transport behavior and leads to the development of a pronounced Gaussian-like density peak (with density excess of a few compared to its surroundings) in the neutrino-cooling layer. This local density maximum has a strong influence on the integrated energy loss by neutrino emission. With better zoning the peak becomes better resolved and even grows in size. Convergence of 1D results seems to be achieved for ⩾600 radial zones, but in two dimensions and three dimensions the artificial density peak prevents numerical convergence for all employed radial zonings (cf. Section 5). Figure 2 also displays some data points for multi-dimensional models with better angular resolution (all of them, however, computed with 400 radial mesh points). These will be discussed in Section 5.

A more detailed analysis, which we will report on below and in the Appendix, reveals that the exact values of the critical luminosities as well as the detailed slope of the critical curves seem to depend strongly on the employed description of neutrino effects, whose implementation is subject to a significant degree of arbitrariness if detailed neutrino transport is not included in the modeling (cf. the discussion in Section 2). The fact that Murphy & Burrows (2008) found fairly good overall agreement between their critical relations $L_\nu (\dot{M})$ and those obtained by Burrows & Goshy (1993) is likely to be linked to a basically similar treatment of the neutrino effects.

Before ∼0.2 s after bounce the mass accretion rate in the case of the 15 M_☉ progenitor changes much more rapidly than during the subsequent evolution (Figure 1). For the corresponding $\dot{M}$ values in excess of ∼0.31 M_☉ s⁻¹, the accretion shock is therefore not as close to steady-state conditions as later on. We see a steep rise of our critical curves at $\dot{M} \gtrsim 0.31\,M_\odot$ s⁻¹ (t_exp ≲ 0.25 s), which is a very prominent difference compared to the results of Burrows & Goshy (1993), who assumed steady-state accretion, but in particular also compared to the hydrodynamic results of Murphy & Burrows (2008) and Nordhaus et al. (2010) even in the 1D case. In order to explore possible reasons for this difference, we performed 1D simulations in which the neutrino treatment is copied from Nordhaus et al. (2010) as closely as possible (i.e., the reduction factor of 2.7 in the exponent of $e^{-\tau _{\mathrm{eff}}}$ is not applied and deleptonization is taken into account by using a Y_e(ρ) relation, but not the corresponding entropy changes proposed by Liebendörfer 2005; for more details on these results, see the Appendix). These runs reveal that the steep rise of our $L_\nu (\dot{M})$ curves is caused by a strong increase of the neutrino-cooling rate with higher values of $\dot{M}$, in particular when we apply our neutrino treatment. The corresponding energy losses inhibit explosions for low values of the driving luminosity. The stronger cooling is linked mainly to our reduction of the effective optical depth τ_eff, which we had to apply in order to reconcile the mass accretion rates and explosion times with the lowest driving luminosities for which Murphy & Burrows (2008) had obtained explosions for the 11.2 and 15 M_☉ stars (cf. Section 2). For example, in the case of the 15 M_☉ progenitor a driving luminosity of $L_{\nu _e} = 3.1\times 10^{52}$ erg s⁻¹ triggers an explosion at t_exp ≈ 250 ms per baryon (p.b.) and $\dot{M}_{\mathrm{exp}}\approx 0.32\,M_\odot$ s⁻¹ (Table 1 and Figure 2), whereas with an implementation of neutrino effects closer to that of Nordhaus et al. (2010) the explosion sets in at t_exp ≈ 120 ms p.b. and $\dot{M}_{\mathrm{exp}}\approx 0.8\,M_\odot$ s⁻¹ (see Figure 19 in the Appendix). Shortly before this moment (at 75 ms after bounce) the total energy loss by neutrino cooling is about 10 times lower with the scheme of Nordhaus et al. (2010) than with our neutrino implementation. The latter yields an integrated energy-loss rate of ∼9 × 10⁵² erg s⁻¹ and significant cooling even at densities between 10¹² and 10¹³ g cm⁻³, where the Nordhaus et al. (2010) treatment shows essentially no cooling. Neither the magnitude of the total neutrino-energy loss rate nor the region of energy extraction with our modeling approach are implausible and in disagreement with detailed transport simulations during a stage when the mass accretion rate still exceeds 1 M_☉ s⁻¹ (cf., e.g., Figure 20 in Buras et al. 2006b). In contrast, the Nordhaus et al. (2010) treatment appears to massively underestimate the neutrino-energy extraction from the accretion flow during this time.

These findings demonstrate that the results of the critical $L_{\nu _e}(\dot{M})$ relation in one dimension can be quantitatively as well as qualitatively different with different approximations of neutrino heating and in particular of neutrino cooling. Moreover, this gives reason for concern that the differences of the critical explosion conditions for 2D and 3D simulations seen by Nordhaus et al. (2010) might have been connected to their treatment of the neutrino physics, in particular also because the decrease of the critical luminosity from 2D to 3D they found was only 15%–25%, which is a relatively modest change (much smaller than the 1D–2D difference) and thus could easily be overruled by other effects. Our results for 2D and 3D simulations with a different implementation of neutrino sources confirm this concern.

In the Nordhaus et al. (2010) paper the average entropy of the matter in the gain region, 〈s(t)〉, was considered to be a suitable diagnostic quantity that reflects the crucial differences of 1D, 2D, and 3D simulations concerning their relevance for the supernova dynamics. While in the spherically symmetric case accreted matter moves through the gain region on the shortest, radial paths, non-radial motions can increase the time that shock-accreted plasma can stay in the gain layer and absorb energy from neutrinos. This can raise the mean entropy, internal energy, and pressure in the postshock region and thus support the revival of the stalled supernova shock. This seems to happen in the simulations by Nordhaus et al. (2010), who found that turbulent mass motions in the 3D case can even improve the conditions for the neutrino-heating mechanism compared to the 2D case. A crude interpretation of this difference can be given by means of random-walk arguments, considering the mass motions in convective and turbulent structures as diffusive process in the gain layer (Murphy & Burrows 2008). The question, however, is whether this effect is a robust 2D–3D difference and whether it is the crucial key to successful explosions by the neutrino-heating mechanism.

Our results at least raise doubts. Figure 4 displays the time evolution of the mean entropy in the gain layer for 11.2 and a 15 M_☉ models computed with driving luminosities near the minimum value for which we obtained explosions. While the 11.2 M_☉ model explodes with a luminosity of 1.0 × 10⁵² erg s⁻¹ for all tested resolutions in more than one dimension despite minimal differences between the values of 〈s(t)〉 compared to the 1D counterpart, the 15 M_☉ progenitor develops an explosion for the chosen luminosity of 2.1 × 10⁵² erg s⁻¹ only in the case of higher-resolution 2D runs (this will be further discussed in Section 5). These successful cases, however, do not stick out by especially high values of 〈s(t)〉. On the contrary, they even have lower mean entropies than the unsuccessful 3D models! It is obvious that Figure 4 does not exhibit the clear 1D–2D–3D hierarchy visible in Figure 5 of the Nordhaus et al. (2010) paper, in which it was found to closely correlate with the explosion behavior of their models. Instead, the differences between simulations in the different dimensions are fairly small, and even 2D flows, which unquestionably allow for explosions also when none happen in one dimension, do not appear more promising than the 1D case in terms of the average entropy of the matter in the gain layer.² Similarly, 3D models possess slightly (insignificantly?) higher values of 〈s(t)〉 but do not show a clear tendency of easier explosions, in particular not the better resolved simulations (see Section 5).

The radial entropy structures of 1D, 2D, and 3D runs for both progenitors, shown in Figures  5 and 6 once before an explosion begins and at another time around the onset of an explosion in at least one of the runs, demonstrate that low-entropy accretion downdrafts and high-entropy rising plumes of neutrino-heated plasma lead to large local variations of the entropy per baryon of the matter in the gain layer (scatter regions in the plots). However, the mass-weighted angular averages of the entropies reveal much smaller differences between the 1D and 2D cases than are visible in Figure 4 of the Nordhaus et al. (2010) paper and in Figure 13 of Murphy & Burrows (2008), and exhibit no obvious signs of more advantageous explosion conditions in the 3D cases compared to the 2D case. The noticeable differences in the radial profiles seem to be insufficient to cause major differences in the mean entropies computed by additional radial averaging (see Figure 4).

How can this discrepancy compared to Nordhaus et al. (2010) and Murphy & Burrows (2008) be explained, and how can one understand the fact that 2D effects play a supportive role for neutrino-driven explosions? We will return to these questions in Section 6, but before that we shall present our results of multi-dimensional simulations with varied resolution in the following section.

In order to test whether our results for the multi-dimensional models depend on the agreeably moderate 3° angular resolution used in the standard runs, we performed a large set of simulations with finer grid spacing especially in the angular directions, but also in radial direction. For this purpose we concentrated on cases around the minimum values of the driving luminosity that triggered explosions of both progenitors in our standard runs. The results are listed in Table 2. They indicate a very interesting trend: 2D models with finer angular zoning tend to explode more readily, whereas better angular resolution in 3D simulations turns out to have the opposite effect.

In the case of the 11.2 M_☉ progenitor, for example, the 3D explosion found to set in at about 730 ms p.b. with an angular zone size of 3° for $L_{\nu _e} = 0.9\times 10^{52}$ erg s⁻¹ cannot be reproduced with an angle binning of 2°. Moreover, a luminosity of $L_{\nu _e} = 1.0\times 10^{52}$ erg s⁻¹ leads to an explosion of the 3D model at ∼540 ms after bounce with a 3° grid, but ∼35 ms later when 2° are used. The corresponding 2D models show the inverse trend as visible in the top panel of Figure 7.

Note that the average shock radii plotted in Figure 7 as well as Figure 3 exhibit alternating periods of increase and decrease in particular in 2D simulations. Such features are a consequence of the strong sloshing motions of the accretion shock and of the associated time-dependent, large global shock deformations, which are typical of violent activity by the SASI. In 3D, the corresponding wiggles and local maxima of the average shock trajectory are much less pronounced. A measure of the degree of shock asphericity, irrespective of the relative weights of different spherical harmonics components, is the standard deviation of the shock radius defined by $\sigma _\mathrm{S} \equiv \sqrt{\vphantom{A^A}\smash{\hbox{${({1}/{4\pi })\oint \mathrm{d}\Omega \, [R_\mathrm{S}(\bm{\Omega })-\langle R_\mathrm{S}\rangle ]^2}$}}}$. The standard deviations corresponding to the average shock radii of Figure 7 are plotted in Figure 8, which confirms the mentioned difference between 2D and 3D runs.

Zoom In Zoom Out Reset image size

In spite of this 2D–3D difference of the shock asphericity, an inspection of cross-sectional snapshots of our best-resolved simulations of the 11.2 M_☉ progenitor with an electron–neutrino luminosity of $L_{\nu _e} = 1.0\times 10^{52}$ erg s⁻¹ reveals that the sizes of the convective plumes and the structure of the neutrino-heated postshock layer are fairly similar in the 2D and 3D cases before explosion (which in both models develops shortly after 500 ms). In Figure 9, it is difficult to judge by eye inspection whether the displayed simulation was conducted in two dimensions (left half-panels) or three dimensions (right half-panels). Even after the explosions have taken off, the global deformation of the shock in both cases is not fundamentally different in the sense that low-order spherical-harmonics modes (dipolar and quadrupolar components) determine the global asymmetry of the shock surface and in particular of the distribution of downdrafts and expanding bubbles in the gain region (see Figure 9, lower right panel, and Figure 10). At a closer inspection one can note some secondary differences in the morphology of the convective and downflow features. Despite the same angular resolution the images of Figures 9 and 10 reveal more small structures in the 3D case compared to the 2D data, which appear more coherent, smoother, and less fragmented into finer substructures and filaments. We will refer to this observation in Section 6.

Our 15 M_☉ runs with varied resolution confirm the trends seen for the 11.2 M_☉ progenitor. For a neutrino luminosity of $L_{\nu _e} = 2.1\times 10^{52}$ erg s⁻¹, for which neither 2D nor 3D simulations with standard resolution produce an explosion, we find that 2D models with angular binning of 15 or better do explode, whereas explosions in 3D models cannot be obtained with angular zones in the range from 15 to 3° (Figure 7, middle panel; Table 2). The 2D simulations exhibit violent SASI sloshing motions and the quasi-periodic appearance of large shock asymmetries (Figure 8, middle panel), and the 2D model with 15 angular zoning explodes with a huge prolate deformation (Figure 11). A similar behavior is seen for the 15 M_☉ runs with $L_{\nu _e} = 2.2\times 10^{52}$ erg s⁻¹: while all 2D models computed with angular zone sizes between 05 and 3° explode, we observe an explosion for the 3D calculation with 3° but none for the case with 2° angular binning (Table 2 and Figure 7, bottom panel). It is highly interesting that the 3° case, where the explosion occurs more readily in three dimensions than in two dimensions, is associated with a large asphericity of the supernova shock at the time the 3D run begins to develop the successful blast (Figure 8, bottom panel). Note again that the structures of the higher-resolved 3D model in Figure 11 reveal finer details and fragmentation into smaller filaments than the corresponding 2D simulation, despite both having the same zone sizes in the angular directions.

The data listed in Table 2 contain the clear message that 2D models with better angular resolution usually develop explosions earlier in contrast to 3D runs, which explode later or not at all when the angular zoning is finer. There can be 2D exceptions to the general trend (e.g., the 15 M_☉ cases with $L_{\nu _e} = 2.1\times 10^{52}$ erg s⁻¹ and 05 and 1° resolution for 400 radial zones; see also Figure 7, middle panel), which are either affected by the difficulty to exactly determine the onset of the blast in cases with a highly deformed shock, or which can be stochastic outliers associated with the chaotic processes leading to the explosion. Note in this context that at late times $\dot{M}(t)$ is very flat and therefore differences in t_exp correspond to only small differences in $\dot{M}_\mathrm{exp}$. For this reason the beginning of the explosion can be shifted by minor perturbations, e.g., connected to stochastic fluctuations. It is also possible that for special circumstances the symmetry axis of the 2D geometry has an influence on such a non-standard behavior because of its effect to redirect converging flows outward or inward and thus to have a positive feedback on the violence of the SASI activity.

Improved radial resolution for fixed angular grid turns out to have a negative influence on the possibility of an explosion also in multi-dimensional simulations. 2D runs with better radial zoning (600 or even 800 instead of 400 radial zones) fail to develop explosions or explode significantly later than their less well-resolved counterparts (see the 15 M_☉ results for $L_{\nu _e} = 2.1\times 10^{52}$ erg s⁻¹ and 2.2 × 10⁵² erg s⁻¹ in Table 2). In general, in 1D, 2D, and 3D runs, improved radial resolution shifts the onset of the explosion to later times monotonically. The only successful 3D model in the set that is useful for the present discussion, a 15 M_☉ run with $L_{\nu _e} = 2.2\times 10^{52}$ erg s⁻¹ and 3° angular resolution, supports our findings in two dimensions. It shows an increasingly delayed explosion for better radial resolution, although the results with 600 and 800 radial zones appear to be nearly converged. All nonexploding 3D models do not yield successes also with higher radial resolution.

The conclusion that good radial resolution is very important for reliable results, in particular when the explosion is "marginal," would not be a surprise because Sato et al. (2009) have pointed out the importance of the radial zoning close to the neutron star and around the supernova shock in order to accurately capture the entropy and vorticity production at the shock and to determine growth times and oscillation frequencies of the SASI. The latter is unquestionably an essential ingredient for the success of the neutrino-driven mechanism in our 2D runs and it may as well be a crucial component for the mechanism to work in 3D runs. As mentioned in Section 4, however, the sensitive influence of the radial zoning in the discussed model set is mainly a consequence of an artificial density peak developing in the neutrino-loss region because of the use of the simplified neutrino-cooling treatment. The narrow shape of this numerical artifact in the density structure, which enhances the energy emission by neutrinos, can even cause a dependence of the results on the particular choice of the grid-cell locations. Different from the 1D runs the results of multi-dimensional simulations with more than ∼600 radial zones do not seem to converge. Since the local density maximum lies between 10¹² g cm⁻³ and 10¹³ g cm⁻³ in the core that is treated spherically symmetrically in our simulations, we interpret this phenomenon as the consequence of a subtle feedback between higher zoning and cooling strength on the one hand and multi-dimensional processes in the accretion layer on the other hand. Because of the artificial nature of the underlying density feature, however, we did not further explore this finding.

Finally, we remark that prior to our present work, Scheck (2007) already performed resolution studies with a large set of 2D simulations, in which he varied the lateral zone width between 05 and 4°. In addition, he conducted three 3D simulations with angular bin sizes of 2°–4°. However, instead of the highly simplified heating and cooling description used by us, he employed the much more sophisticated approximation for gray neutrino transport described in detail in Scheck et al. (2006). This approximation included, e.g., the feedback of accretion on the neutrino emission properties and on the corresponding energy and lepton-number transport by neutrinos and antineutrinos of all flavors, as well as a more elaborate description of neutrino-matter interactions in detailed dependence on the thermodynamical state of the stellar plasma. Scheck (2007) was not interested in a systematic exploration of the critical explosion condition, but his project was focused on investigating the possibility of hydrodynamic pulsar kicks by successful asymmetric explosions. Despite the grave differences of the neutrino treatments and numerical setups, the results obtained by Scheck (2007) are compatible with our present findings: 2D simulations with higher resolution turned out to yield explosions significantly earlier and thus also more energetically than the low-resolution runs. Within the tested range of angular resolutions, Scheck (2007) did not observe any significant differences between his 3D models. This, however, may just be a consequence of the fact that the models were clearly above the threshold conditions for an explosion and did not linger along the borderline between blast and failure.

In this section, we discuss the meaning of our results in comparison to previous studies and present an interpretation that could explain the main trends found in our multi-dimensional simulations with varied resolution.

6.1. Variation with Dimension

In Sections 4 and 5, we have reported that our simulations do not support the central finding by Nordhaus et al. (2010) that the tendency to explode is a monotonically increasing function of dimension. While we confirm more favorable explosion conditions in two dimensions than in one dimension we do not observe that in three dimensions considerably lower driving luminosities are needed for a success of the neutrino-driven mechanism than in two dimensions. Moreover, we cannot confirm the finding by Nordhaus et al. (2010) that the mass-weighted average of the entropy per nucleon in the gain region, 〈s(t)〉, is a quantity that is suitable as an indicator of the proximity of models to an explosion and thus can serve as an explanation of differences between 1D, 2D, and 3D simulations. In particular, our 3D models turned out to have slightly higher mean entropies than corresponding 2D cases (Figure 4) without developing better explosion conditions. This raises the question of why our models, and multi-dimensional simulations in general, have produced successful explosions by the neutrino-heating mechanism when corresponding 1D models fail.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is by no means obvious that 〈s(t)〉 should increase in the gain layer in the multi-dimensional case. While neutrino-energy deposition naturally leads to a rise of the entropy of the heated gas, the averaging process over the volume of the gain layer also encompasses the downdrafts carrying cool matter from the postshock region toward the gain radius and the neutron star. These downdrafts are much denser, they are hardly heated by neutrinos because of their extremely rapid infall, and they can contain more mass than the surrounding, dilute bubbles that are inflated by the expanding, neutrino-heated plasma. It is therefore not clear that the spatial (mass-weighted) average 〈s(t)〉 grows in multi-dimensions compared to 1D runs.

Moreover, it is not even clear that convective overturn in the gain layer must lead to an average entropy of the neutrino-heated gas itself that is higher than in 1D simulations. Different from the 1D case high-entropy matter becomes buoyant and begins to float in the multi-dimensional environment. Thus, the heated gas is quickly carried away from the vicinity of the gain radius, where neutrino-energy deposition is maximal, to larger radii. Such dynamics of the gas can limit well the amount of energy and entropy that is stored in individual chunks of matter. Little, if any of the gas is subject to multiple overturn cycles bringing the gas close to the gain radius more than once as suggested by the "convective engine" picture of Herant et al. (1994) but questioned by Burrows et al. (1995). Instead, the majority of the heated gas either expands to larger distances, pushing shock expansion, or, in the disadvantageous case, is swept back below the gain radius (e.g., by large-amplitude sloshing motions of the shock), where it loses its energy again by efficiently reradiating neutrinos.³

Our results imply that the dominant effect that makes the multi-dimensional case more favorable for an explosion than spherical symmetry is associated with an inflation of the shock radius and postshock layer, driven by the buoyant rise and expansion of the plumes of neutrino-heated plasma. In the course of the postshock volume becoming more extended, the integrated mass M_gain in the heating layer increases compared to the 1D case. This can be seen in Figure 12, which displays the mass in the gain layer as function of postbounce time for the 11.2 and 15 M_☉ runs with the different neutrino luminosities and resolutions already shown in Figures 7 and 8. While the mass-averaged entropy 〈s〉 in the gain region hardly changes, the integral value of the entropy, M_gain〈s〉, clearly increases with the models coming closer to explosion. This dependence is particularly well visible when 2D and 3D models with different resolutions are compared with each other.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The longer dwell times of matter in the gain layer of 2D simulations observed by Murphy & Burrows (2008), which correspond to the advection times τ_adv evaluated by Buras et al. (2006a) and Marek & Janka (2009) (though different ways of calculation have been considered, in particular for the multi-dimensional case), are a manifestation that a growing mass accumulates in the postshock region to get energy loaded by neutrino absorption and to finally drive the successful supernova blast. In near-steady-state conditions the mass accretion rate through the gain layer is equal to the mass infall rate $\dot{M}$ ahead of the shock, where it is determined by the core structure of the progenitor star. Since $\tau _\mathrm{adv} \approx M_\mathrm{gain}/\dot{M}$ (cf. Marek & Janka 2009), a larger value of τ_adv correlates with a higher mass M_gain. Accordingly, the total net heating rate Q_gain and thus the heating efficiency $\epsilon \equiv Q_\mathrm{gain}/(L_{\nu _e}+L_{\bar{\nu }_e}) = Q_\mathrm{gain}/(2L_\nu)$ of the gas residing in the gain layer is also higher for models that develop an explosion (see Figure 13 and Murphy & Burrows 2008). Figure 14 shows the ratio of the advection (dwell) timescale to the neutrino-heating timescale in the gain layer (our evaluation employs the Newtonian analog of the formulas given in Müller et al. 2012) for the models also displayed in Figures 12 and 13. The same trends as in the previous images can be seen. Models closer to an explosion exhibit higher values of the timescale ratio. The ratio approaches unity, indicating the proximity to a runaway instability, roughly around the time when we define the onset of an explosion (i.e., when the average shock radius of Figure 7 passes 400 km). The quality of the coincidence of these moments, however, differs from model to model and depends on the exact definition used for the timescales and the accuracy at which the relevant quantities can be evaluated in highly perturbed flows. Nevertheless, the evolutionary trend of the timescale ratio (increasing or decreasing) is suggestive for whether a simulation run leads to a successful explosion or failure.

A larger mass in the gain layer and higher total net energy deposition rate are therefore better indicators of the proximity of our models to explosion than the mean entropy of the gas in this region, which does not exhibit the 1D–2D–3D hierarchy with dimension found previously by Nordhaus et al. (2010). As discussed in Section 4 the main reason for this discrepancy is most probably the different treatments of neutrino lepton-number losses and our consequential recalibration of the energy source terms. This leads to significantly higher energy drain from the cooling layer in our simulations. While this hypothesis is supported by 1D tests that we conducted, we cannot be absolutely certain that no other effects play a role for the discrepancies between our results and those of Nordhaus et al. (2010), because detailed cross-comparisons are not available and our knowledge of the details of the implementation of neutrino effects by Nordhaus et al. (2010) may be incomplete. Other potential reasons for differences may be connected to the hydrodynamics scheme (Prometheus with a higher-order Godunov solver and directional splitting versus Castro with unsplit methodology; Almgren et al. 2010), the employed grid (polar coordinates versus structured grid with adaptive refinement by a nested hierarchy of rectangular grids), potentially—though not very likely—the use of a 1D core above 10¹² g cm⁻³ in our simulations, or differences in the exact structure and properties of the infall region upstream of the stalled shock as a consequence of different treatments of the collapse phase until 15 ms after core bounce (due to full neutrino transport plus a nuclear-burning approximation in the Prometheus-Vertex code versus the simple deleptonization scheme of Liebendörfer 2005) or of different seeding of non-radial hydrodynamic instabilities (in our case by imposed, small random seed perturbations of the density), or linked to differences of the low-density EoS outside of the application regime of the Shen et al. (1998) EoS.

Despite these uncertainties about the exact cause of the differences, whose ultimate elimination will require systematic and time-consuming studies, our results, as they are, send a clear message. The outcome of the 1D–2D–3D comparison and the effects of the third dimension advertised by Nordhaus et al. (2010) "as a key to the neutrino mechanism of core-collapse supernova explosions" are not at all robust results. Instead, the exact slope of the critical explosion condition $L_\nu (\dot{M})$, its location, and its shift with dimension, as well as the existence of a 1D–2D–3D hierarchy of the mass-averaged entropy in the gain layer seem to depend sensitively on subtle details of the neutrino treatment or other numerical aspects of the simulations.

6.2. Resolution Dependence

Let us now turn to the second, highly interesting question in connection with our set of models, namely, to the resolution dependence of our results. Our set of simulations performed with different angular binnings reveals that quantities that turned out to diagnose healthy conditions for an explosion, i.e., the growth of the average shock radius, the degree of shock deformation, or the mass and total heating in the gain layer (but not the mass-averaged entropy of the matter in the gain region), show a clear dependence on the angular zoning (see Figures 7, 8, 12, and 13 in contrast to Figure 4). In particular, 2D models with better angular resolution exhibit more favorable conditions and explode more readily (in agreement with results obtained by Scheck 2007 with a more sophisticated treatment of neutrino transport than the simple heating and cooling source terms applied in our investigation), whereas 3D models obey the opposite behavior. What does the resolution dependence of our simulations tell us about the mechanism leading to explosions in our models? And how can we understand the puzzling finding that 2D and 3D runs follow opposite trends when the angular resolution is refined?

We interpret this as a manifestation of two aspects or facts.

• 1.  
  The success of our models, at least in the neighborhood of the explosion threshold, is fostered mainly by large-scale mass flows as associated with strong SASI activity, but not by enhanced fragmentation of structures and vortex motions on small spatial scales.
• 2.  
  Our resolution study reflects the consequences of the turbulent energy cascade, which redistributes energy fed into the flow by external sources in opposite directions in two dimensions and three dimensions while in three dimensions the turbulent energy flow goes from large to small scales, it pumps energy from small to large spatial scales in two dimensions.

Point (1) is supported by the kinetic energies of non-radial mass motions in the gain layer of the 2D and 3D models plotted in Figure 15. From this picture it is obvious that in the case of successful models the angular kinetic energy is higher and shows an overall trend of growth in time until the blast has taken off. Moreover, the spiky maxima and minima of quasi-periodic variations, which are indicative of the presence of low-order SASI modes, are significantly larger for exploding models. This not only holds for 2D models, whose lateral kinetic energies exhibit variation amplitudes of several 10% and partially up to even ∼50% of the time-averaged value, but is also true for 3D models, although in this case the amplitudes are generally smaller and the non-radial kinetic energy is split essentially equally into lateral and azimuthal contributions. When comparing successful 2D runs with those in 3D runs, our studies suggest that in both cases the shock exhibits a growing degree of asphericity (expressed by the standard deviation of the shock deformation plotted in Figure 8) when the explosion is approached, and the kinetic energy of non-radial mass motions reaches roughly the same magnitude (Figure 15), at least for models near the explosion threshold.

Actually, a variety of observations can be interpreted as support of the hypothesis that flows on the largest possible scales rather than on small scales play a crucial role for the success of the neutrino-heating mechanism in our simulations.

• 1.  
  The strength of low-mode SASI activity in 2D models as indicated by growing fluctuations of the angular kinetic energy and of the shock deformation (Figures 15 and 8) increases with higher resolution in clear correlation with an earlier onset of the explosion (Figure 7). Stronger SASI activity obviously facilitates explosions, which is visible by a growing average shock radius as well as larger mass and higher total heating rate in the gain layer.
• 2.  
  Exploding 2D models as well as 3D models exhibit large shock deformation at the time of explosion (although the relative asphericity σ_S/R_S of the shock surface is somewhat smaller in 3D models than in 2D models; see Figures 7 and 8).
• 3.  
  More fine structure on small spatial scales, which can be seen in 3D models computed with higher resolution in Figures 9–11, does not imply improved conditions for an explosion.
• 4.  
  Exploding 2D models are not connected with the highest mean entropies in the gain region (Figure 4). This fuels doubts in a random-walk picture where turbulent vortex motions on small scales enlarge the residence time of matter in the gain layer (Murphy & Burrows 2008) and thus could allow for more energy absorption of such mass elements from neutrinos. If this effect occurred, it does not concern a major fraction of the mass in the gain region.

Point (2) is the only plausible argument we can give for explaining the opposite response to higher angular resolution that we discovered in our 2D and 3D simulations. The sequence of 2D runs with gradually reduced lateral zone sizes reflects the growing violence of large-scale flows by higher fluctuation amplitudes of the kinetic energy in the gain layer (Figure 15) and larger temporal variations of the average shock radius and shock deformation (Figures 7 and 8). In contrast, more energy on small spatial scales in the 3D case manifests itself by a progressing fragmentation of the flow, leading to a growing richness of vortex structures and finer filaments in the case of 3D models with smaller angle bins (Figures 9–11). As a consequence, 3D models with higher angular resolution become more similar to the 1D case in various quantities that we considered as explosion indicators, see, e.g., Figures 7, 8, 12, and 13.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Both the powerful coherent mass motions of the SASI layer in 2D and the vivid activity in small vortex structures in the 3D environment are fed by two external sources which supply the postshock layer with an inflow of fresh energy: (1) gravitational potential energy that is released by the continuous stream of matter falling through the accretion shock and (2) energy deposition by neutrinos. The energy stored in the fluid is then redistributed toward small or large scales according to the turbulent cascades characteristic of 2D and 3D environments.

Direct evidence for the action of different turbulent energy cascades in 2D and 3D can be obtained by considering the energy spectrum E(k) of turbulent motions as a function of wavenumber k in the gain region. The spectral shape of E(k) can already be adequately established by considering only the azimuthal velocity v_θ at a given radius using a decomposition into spherical harmonics Y_lm(θ, ϕ):

$$\begin{equation} E(l) = \sum ^{l}_{m=-l}\left|\int _{\Omega }Y_{lm}^{*}(\theta ,\phi) \, \sqrt{\rho } \, v_{\theta }(r,\theta ,\phi) \, d\Omega \right|^{2} \,. \end{equation} \tag{ 8 }$$

Here, the velocity fluctuations have been expressed in terms of the multipole order l instead of the wavenumber k. A summation over the energies of modes with the same l has been carried out,⁴ and in order to obtain smoother spectra, we average E(l) over 30 km in radius and over 10 time steps. One expects that the resulting spectrum E(l) directly reflects the properties of E(k) such as the slopes in different regimes of the turbulent cascade.⁵ The computed spectra E(l) (Figure 16) indeed confirm the predictions from 3D and (planar) 2D turbulence theory, at least for sufficiently high multipole order l.

Zoom In Zoom Out Reset image size

In 3D (right panel), a power-law spectrum with E(l)∝l^−5/3 (Landau & Lifshitz 1959) develops at intermediate wavenumbers as the resolution is increased, reflecting the transfer of energy from large to small scales in a forward cascade until dissipation takes over at large l. At high resolution, the energy contained in small-scale disturbances increases, as the dissipation range moves to larger l. One observes that the 5/3 power law is broken at low l (l ≲ 10), suggesting that kinetic energy is injected into the flow at wavenumbers l ≈ 10, i.e., at scales typical for growing convective plumes.

By contrast, the power-law dependence E(l)∝l^−5/3 approximately holds for l ≲ 10 in 2D as a result of the reverse energy cascade (Kraichnan 1967). The energy injected at l ≈ 10 is therefore not transferred to the dissipative range; only enstrophy (the squared vorticity of the velocity field) is transported in a forward cascade with a different power-law index (E(l)∝l⁻³).

This appears to be a natural explanation for the predominance of large-scale and small-scale structures in 2D and 3D models, respectively. Moreover, this picture suggests that as dissipation affects the "injection scale" at l ≈ 10 less with increasing resolution, the differences between 2D and 3D models become more pronounced with finer grid zoning.

Our analysis of the spectral properties of turbulence thus further strengthens our view that the trends seen in our simulations strongly suggest that non-radial kinetic energy available on large scales, not on small scales, assists the development of an explosion by the neutrino-heating mechanism. This explains why 2D models with higher angular resolution tend to explode earlier and thus at higher values of the mass-accretion rate than less resolved models. On the other hand, the energy "drain" by vortex motions on ever smaller scales—with the same reservoir of pumping energy per unit mass being available from accretion and neutrino heating—disfavors explosions in better resolved 3D models.

We therefore conclude that the key to the mechanism of core-collapse supernova explosions seems intrinsically and tightly linked to the question how much kinetic energy of the matter in the gain region can be accumulated in non-radial fluid motions on the largest possible scales, i.e., in the lowest-order spherical harmonics modes of non-radial hydrodynamic instabilities. The predominant growth of such flows is typical of SASI activity, whose lowest-order spherical harmonics modes possess the highest growth rates (Blondin et al. 2003; Blondin & Mezzacappa 2006b; Foglizzo et al. 2006, 2007; Ohnishi et al. 2006). Strong SASI motions drive shock expansion, increase the gain layer and its mass content, allow a larger fraction of the accreted matter to stay in the gain layer and be exposed to efficient neutrino heating, and thus aid the development of an explosion (Scheck et al. 2008; Marek & Janka 2009). However, our models do not show a systematic trend of higher average entropies of the matter in the gain layer for models closer to explosion. Instead, we find that such models have larger mass, larger non-radial kinetic energy, larger total neutrino-heating rate, and larger total entropy in the gain layer.

We have performed a systematic study of the postbounce evolution of supernova cores of 11.2 and 15 M_☉ and their explosion by the neutrino-heating mechanism in one, two, and three dimensions, employing simple neutrino-cooling and heating terms with varied values of the driving luminosity. We conceptually followed previous studies by Murphy & Burrows (2008) and Nordhaus et al. (2010), but did not apply the deleptonization treatment that they adopted from Liebendörfer (2005), who introduced it for an approximative description of neutrino losses during the infall phase until core bounce. We argued (Section 2) that this approximation—with or without the source term proposed by Liebendörfer (2005) to account for entropy generation in neutrino-electron scatterings—does not provide a suitable treatment of the evolution of the electron abundance after core bounce. Therefore, we did not consider changes of the net electron fraction Y_e of the stellar plasma at times later than 15 ms after bounce, up to which the collapse was followed with the Prometheus-Vertex code including full neutrino transport. While ignoring Y_e changes subsequently is certainly not a good approximation, it is not necessarily more unrealistic than describing the lepton-number evolution during the accretion phase of the stalled shock by the scheme of Liebendörfer (2005). As a consequence, we had to replace an exponential factor $e^{-\tau _{\mathrm{eff}}}$, which was introduced in an ad hoc way by Murphy et al. (2009) and Nordhaus et al. (2010) to damp the neutrino source terms at high optical depths τ_eff, by $e^{-\tau _{\mathrm{eff}}/2.7}$ in order to reproduce the minimum luminosity found to yield explosions in the 1D simulations by Murphy & Burrows (2008) and Nordhaus et al. (2010). This modification led to enhanced neutrino losses in the cooling layer, which were better compatible with total energy loss rates found in simulations with detailed neutrino transport, e.g., in Buras et al. (2006b) and is responsible for some of the findings and differences discussed in Section 4.

Our results and conclusions can be briefly summarized as follows.

• 1.  
  We cannot reproduce the exact slopes and relative locations of the critical curves $L_\nu (\dot{M})$ of 1D, 2D, and 3D simulations found by Nordhaus et al. (2010). While our results confirm the well-known fact that explosions in two dimensions occur for a lower driving luminosity L_ν than in one dimension when the mass accretion rate $\dot{M}$ is fixed, we cannot discover any significant further reduction when we go from 2D to 3D models.
• 2.  
  We cannot confirm that the mass-averaged entropy of the matter in the gain region, 〈s〉, is a good diagnostic quantity for the proximity to an explosion. As we argued in Section 6.1, it is neither clear nor necessary that 〈s〉 is higher for cases where explosions are obtained more readily. Our successful 2D models do not exhibit larger mean entropies than the corresponding 1D cases, which fail to explode. Instead, we observed that the total mass, total entropy, total neutrino-heating rate, and the non-radial kinetic energy in the gain layer are higher in cases that develop an explosion.
• 3.  
  We conclude that the tendency for an explosion as a monotonically increasing function of dimension as well as the 1D–2D–3D hierarchy of 〈s(t)〉 found by Nordhaus et al. (2010) are not robust results. They seem to be sensitive to subtle differences of the approximations of neutrino effects (and/or to other differences in the numerical treatments of the models). It is therefore unclear how far studies with radical simplifications of the neutrino physics (without detailed energy and lepton-number source terms and transport; no feedback between accretion and neutrino properties) can yield results that are finally conclusive for the explosion-triggering processes in real supernova cores.
• 4.  
  Increasing the angular resolution we observed a clear tendency of 2D models to explode earlier, in agreement with previous results by Scheck (2007), who employed a more sophisticated treatment of neutrino effects based on the transport approximation described in Scheck et al. (2006). In contrast, 3D models show the opposite trend and in a variety of quantities and aspects become more similar to their 1D counterparts. The easier explosion of the 2D models is connected to an enhanced violence of large-scale mass motions in the postshock region due to SASI activity, whereas 3D models with better angular resolution appear to develop less strength in low-order SASI modes.
• 5.  
  We interpret this finding as a consequence of the opposite turbulent energy cascades in two dimensions and three dimensions. In two dimensions, the energy continuously pumped into the gain layer by neutrino heating and the release of gravitational binding energy flows from small to large scales and thus helps to power coherent mass motions on the largest possible spatial scales. In contrast, in three dimensions this energy seems to instigate flow vorticity and fragmentation of structures on small scales. Evidence for this interpretation is provided by Figure 16.
• 6.  
  We also conclude from our resolution studies that the presence of violent mass motions connected to low-order SASI modes is favorable for an explosion (in agreement with arguments given by Marek & Janka 2009 and Scheck et al. 2008). This is supported by the fact that 2D and 3D models that are closer to explosion show signs of growing power in large-scale mass motions (signaled by growing fluctuations of the kinetic energy of non-radial velocity components) and in particular develop significant shock deformation and global ejecta asymmetries when the explosion sets in.
• 7.  
  We found that higher radial resolution makes explosions more difficult with the setup chosen for the investigated set of models. Higher resolution turned out to prevent explosions or to let them occur later in simulations in 1D, 2D, and 3D. This result could be diagnosed to be a consequence of a local density maximum in the neutrino-cooling layer, which grows with higher resolution and enhances the energy loss by neutrino emission. This density peak, however, is a numerical artifact of the employed simple neutrino-cooling treatment by an analytic source term which is exponentially damped at high densities.

The lack of very precise information on the physics ingredients and their exact implementation, e.g., details of the treatment of neutrino source terms, low-density EoS, and progenitor data when mapped into the simulation and seeded with small random perturbations, as well as a variety of methodical differences like the hydrodynamics scheme, numerical grid, and the use of a 1D core at high densities or not, prevent us from presenting a rigorous proof that could causally link the discrepancies between our results and those of Nordhaus et al. (2010) to one or more well understood reasons. We think that the nagging uncertainties in this context demand a future, involved, collaborative code-comparison project. This will also require considerable amounts of computer time for further 3D simulations, in particular with high resolution, thus needing more computer resources than available to us for the described project.

Despite this deficiency, however, our results suggest that the differences of 3D compared to 2D simulations observed by Nordhaus et al. (2010) are unlikely to be a robust outcome but seem to depend on relevant aspects of the modeling (most probably the neutrino physics, but potentially, and not finally excluded, also technical aspects).

We therefore conclude that the influence of 3D effects on the supernova mechanism is presently not clear. We strongly emphasize, however, that the fact that our results do not corroborate improved explosion conditions in 3D compared to 2D cannot be used as an argument that 3D effects do not facilitate the supernova explosion mechanism or are of minor importance. We just think that in the context of the neutrino-driven mechanism the relevance and exact role of 3D fluid dynamics are not yet understood. We therefore have the opinion that the results obtained by Nordhaus et al. (2010) do not justify their claims that 3D hydrodynamics offers the key to a fundamental understanding of the neutrino mechanism while other physics in the supernova core, like general relativity or the properties of the nuclear EOS, are only of secondary importance. Though this may well be right, such statements at the present time are premature and not supported by solid facts and results.

Our study, however, raises further important questions. How far can our understanding be developed on grounds of modeling approaches that employ radical simplifications of the neutrino physics? Which aspects of the complex interplay between different components of the problem are linked to the essence of how the explosion is triggered by the combination of neutrino-energy supply and non-radial hydrodynamic instabilities? Examples for such mutually dependent components are the neutrino transport and hydrodynamics, the neutron star core evolution and fluid motions around the neutron star, or the mass flux from the accretion shock to the deceleration layer (both being the coupling regions for the advective–acoustic cycle that is thought to be responsible for the SASI growth; e.g., Scheck et al. 2008) and the conditions in the neutrino-heating and cooling layers. Much more work needs to be done to find the answers of these questions.

A major shortcoming of the setup applied in previous works and adopted also in our investigation is the neglect of neutrino cooling and deleptonization inside the proto-neutron star. The employed simple neutrino source terms and their exponential suppression at high optical depths do not allow the neutron star to evolve. Underestimating the neutron star contraction, however, slows down the infall velocities in the postshock layer and thus has disfavorable consequences for the growth of the SASI (Scheck et al. 2008), similar to the effects of reduced neutrino losses in the cooling layer (Scheck et al. 2008) or increased nuclear photodissociation behind the stalled shock (Fernández & Thompson 2009). On the contrary, a slower postshock flow improves the conditions for the growth of convective instability, whose development is supported by a high ratio (larger than ∼3 signals a linearly unstable situation) of the advection timescale through the gain region divided by the inverse of the Brunt–Väisälä frequency as discussed by Foglizzo et al. (2006), Buras et al. (2006a), and Scheck et al. (2008). It is therefore well possible that the weak SASI activity diagnosed in recent 2D and 3D simulations by Burrows et al. (2012)—leading them to the conclusion that dipolar asymmetries are caused by convection rather than the SASI—is an artifact of the approximative treatment of neutrino cooling and the disregard of neutron star contraction in their models. In full-scale supernova models with sophisticated neutrino transport, Buras et al. (2006a) observed differences in the growth conditions for convection compared to the SASI in collapsing stellar cores of a variety of progenitor stars. Also, the high-density EoS and general relativity, influencing the contraction behavior of the nascent neutron star, can make a difference (Marek & Janka 2009; Müller et al. 2012). Claims, based on highly simplified models, that SASI is less important than convection and is at most a minor feature of the supernova dynamics (Burrows et al. 2012) are therefore certainly premature.

Finally, our resolution study suggests that the action of the turbulent cascade in 3D extracts energy from coherent large-scale modes of fluid motion and instead fuels fragmentation and enhanced vortex flows on small spatial scales. At least in our 3D models with better grid zoning the appearance of finer structures in the postshock flow was connected with a tendency of damping the development of explosions. While a finally convincing proof of such a negative feedback may require much better resolved simulations than we presently can afford to conduct (in order to minimize numerical dissipation on small scales), this result implies that good resolution—considerably higher than recently used by Takiwaki et al. (2011), whose 3D simulation had only 32 azimuthal zones (corresponding to a cell size of 1125)—is indispensable to clarify the 3D effects on the explosion mechanism. Moreover, our result points in an interesting direction. Is it possible that the success of the neutrino-driven mechanism in three dimensions is tightly coupled to the presence of violent SASI activity, a connection that was found before—and is confirmed by our present study—to foster explosions in two dimensions? If so, what is the key to instigate such violent SASI motions of the supernova core in three dimensions? Will they occur with a better (more realistic) treatment of the neutrino transport and correspondingly altered conditions in the heating and cooling layers and in the contracting core of the proto-neutron star? Or are they associated with stellar rotation, which even with a slow rate can initiate the faster growth of spiral (nonaxisymmetric) SASI modes (Blondin & Mezzacappa 2006a; Yamasaki & Foglizzo 2008; Iwakami et al. 2009; Fernández 2010)? Or is strong SASI activity in the supernova core triggered by large-scale inhomogeneities in the 3D progenitor star (Arnett & Meakin 2011), which could provide a more efficient seed for SASI growth than the random cell-to-cell small-amplitude perturbations employed in our simulations? Should the presence of large-amplitude SASI mass motions indeed turn out to be the key to the neutrino mechanism in three dimensions, it would mean that neutrino-driven explosions are not only a generically multi-dimensional phenomenon, but one that is generically associated with dominant low-order modes of asymmetry and deformation from the very beginning.

While this paper raises many more questions than it is able to answer, it definitely makes clear that our understanding of the supernova physics in the third dimension is still in its very infancy. A virgin territory with distant horizons lies ahead of us and awaits to be explored.

We are grateful to Lorenz Hüdepohl for his valuable input to different aspects of the reported project and thank Elena Erastova and Markus Rampp (Max-Planck-Rechenzentrum Garching) for their help in the visualization of our 3D data. H.T.J. thanks Rodrigo Fernández, Thierry Foglizzo, Jerome Guilet, and Christian Ott for stimulating and informative discussions. This work was supported by the Deutsche Forschungsgemeinschaft through Sonderforschungsbereich/Transregio 27 "Neutrinos and Beyond," Sonderforschungsbereich/Transregio 7 "Gravitational-Wave Astronomy," and the Cluster of Excellence EXC 153 "Origin and Structure of the Universe." The computations were performed on the Juropa cluster at the John von Neumann Institute for Computing (NIC) in Jülich, partially through a DECI-6 grant of the DEISA initiative, on the IBM p690 at Cineca in Italy through a DECI-5 grant of the DEISA initiative, and on the IBM p690 at the Rechenzentrum Garching.

In this Appendix, we briefly report on our efforts to reproduce the 1D results of Murphy & Burrows (2008) and Nordhaus et al. (2010) for the critical explosion conditions of the 15 M_☉ progenitor, applying a neutrino treatment that was intended to copy the procedure outlined in these publications as closely as possible.

For this purpose we retained the exponential suppression factor $e^{-\tau _\mathrm{eff}}$ of Equations (4) and (5) without a reduction factor of 2.7 in the exponent, and the lepton evolution before and after core bounce was described by employing a predefined Y_e(ρ) relation. We also aimed at reproducing the core infall of the previous works as closely as possible, because the density structure of the infall region ahead of the stalled shock determines the mass-infall rate $\dot{M}(t)$ at the shock, and some differences became visible when we compared our values with those given by Murphy & Burrows (2008) and Nordhaus et al. (2010) (Figure 17). We therefore recomputed the collapse phase from the onset of gravitational instability of the progenitor core through core bounce with the deleptonization scheme of Liebendörfer (2005). Entropy changes were taken into account as suggested by Liebendörfer (2005), but were switched off after core bounce following Murphy & Burrows (2008) and Nordhaus et al. (2010).

Zoom In Zoom Out Reset image size

We tested three different cases for the functional relation Y_e(ρ): first, we used a tabulated result for the Y_e(ρ) evolution as obtained with the Prometheus-Vertex code and state-of-the-art electron-capture rates (Langanke et al. 2003) ("Deleptonization 1"). Second, we applied a Y_e(ρ) table provided by Christian Ott as a co-developer of the CoCoNuT code (http://www.mpa-garching.mpg.de/hydro/COCONUT/). These data are based on collapse simulations with the Vulcan/2D code (Livne et al. 2004) ("Deleptonization 2"). Third, we employed a fitting formula given by Liebendörfer (2005) for the parameters of model G15 ("Deleptonization 3"). All three Y_e(ρ) trajectories are depicted in Figure 18.

Zoom In Zoom Out Reset image size

The three sets of 1D simulations conducted for this Appendix were performed with 800 radial zones. The corresponding critical luminosity curves $L_{\nu _e}(\dot{M})$ are displayed in Figure 19 in comparison to those of Murphy & Burrows (2008) and Nordhaus et al. (2010) and to our results of Figure 2 (for 400 zones, because for this resolution the calibration of the exponential suppression factor for best agreement with the critical curve of Murphy & Burrows 2008 was done). The overall slopes of all three curves are similar but none of them is quantitatively or qualitatively in good agreement with those of Murphy & Burrows (2008) and Nordhaus et al. (2010). Explosions in our simulations occurred significantly more readily (i.e., for lower $L_{\nu _e}$) than in the previous works. This suggests less cooling in our runs, although we made all possible efforts to exactly follow the description of the neutrino treatment in those papers. The steep rise and partly backward bending of our curves for $\dot{M}$ values around 0.2–0.3 M_☉ s⁻¹ can be understood by an inspection of Figure 20, which shows the time evolution of the shock radius for simulations with prescription "Deleptonization 1" for a selection of neutrino luminosity values. One can see that in the case of $L_{\nu _e} = 2.95\times 10^{52}$ erg s⁻¹ the shock makes a larger excursion before it returns again. Its reexpansion, leading to an explosion, therefore happens later than in the model with $L_{\nu _e} = 2.7\times 10^{52}$ erg s⁻¹, where the first shock expansion is much less strong. Correspondingly, the explosion in the former case sets in at a later time and lower mass accretion rate than in the latter case, explaining the backward bending of $L_{\nu _e}(\dot{M})$ in this regime of luminosities and $\dot{M}$. The nearly horizontal parts of the critical curves can be understood by the fact that for such high values of the luminosities the neutrino cooling (with the unmodified $e^{-\tau _\mathrm{eff}}$ suppression factor) is so weak that the explosion sets in very early (see the black line in Figure 20) and therefore for large values of the mass accretion rate. The region between $\dot{M} \approx 0.3\,M_\odot$ s⁻³ and $\dot{M} \approx 0.8\,M_\odot$ s⁻³ is difficult to probe with a stepwise increase of $L_{\nu _e}$, because the mass-accretion rate there changes so rapidly that the shock shows time-dependent dynamics instead of settling into a quasi-steady state.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

None of the critical curves obtained with the direct implementation of the neutrino treatment described in Murphy & Burrows (2008) and Nordhaus et al. (2010) can reproduce the critical luminosity curves reported in these papers reasonably well. We therefore decided to proceed with the modifications described in Section 2.