The Dynamics of Neutrino-Driven Supernova Explosions after Shock Revival in 2D and 3D
MMon. Not. R. Astron. Soc. , 000–000 (0000) Printed July 17, 2015 (MN L A TEX style file v2.2)
The Dynamics of Neutrino-Driven Supernova Explosions afterShock Revival in 2D and 3D
B. M ¨uller (cid:63)
Monash Centre for Astrophysics, School of Physics and Astronomy, 19 Rainforest Walk, Monash University, Victoria 3800, Australia
July 17, 2015
ABSTRACT
We study the growth of the explosion energy after shock revival in neutrino-driven explosionsin two and three dimensions (2D / . M (cid:12) star. The 3D model shows a faster and steadier growth of the explosion energyand already shows signs of subsiding accretion after one second. By contrast, the growth ofthe explosion energy in 2D is unsteady, and accretion lasts for several seconds as confirmedby additional long-time simulations of stars of similar masses. Appreciable explosion energiescan still be reached, albeit at the expense of rather high neutron star masses. In 2D, the bind-ing energy at the gain radius is larger because the strong excitation of downward-propagating g -modes removes energy from the freshly accreted material in the downflows. Consequently,the mass outflow rate is considerably lower in 2D than in 3D. This is only partially com-pensated by additional heating by outward-propagating acoustic waves in 2D. Moreover, themass outflow rate in 2D is reduced because much of the neutrino energy deposition occurs indownflows or bubbles confined by secondary shocks without driving outflows. Episodic con-striction of outflows and vertical mixing of colder shocked material and hot, neutrino-heatedejecta due to Rayleigh-Taylor instability further hamper the growth of the explosion energy in2D. Further simulations will be necessary to determine whether these e ff ects are generic overa wider range of supernova progenitors. Key words: supernovae: general – hydrodynamics – instabilities – neutrinos – radiative trans-fer
After decades of research, the mechanism powering core-collapsesupernovae is still one of the outstanding problems in theoreti-cal astrophysics. The so-called delayed neutrino-driven mechanismcurrently remains the most popular and best explored explanation(see Janka 2012; Burrows 2013 for current reviews) for “ordi-nary” supernova explosions not exceeding ∼ erg. In its mod-ern guise, the neutrino-driven mechanism relies on additional sup-port for shock expansion in the form of hydrodynamic instabilitieslike convection (Herant, Benz & Colgate 1992; Burrows & Fryxell1992; Herant et al. 1994; Burrows, Hayes & Fryxell 1995; Janka& M¨uller 1996; M¨uller & Janka 1997) and the standing accretionshock instability (SASI, Blondin, Mezzacappa & DeMarino 2003;Laming 2007; Foglizzo et al. 2007; Guilet & Foglizzo 2012). In-deed, many successful multi-dimensional simulations of neutrino-driven shock revival (mostly in 2D, i.e. under the assumption ofaxisymmetry) have strengthened our confidence in the neutrino-driven mechanism over the recent years (Buras et al. 2006a; Marek& Janka 2009; Yakunin et al. 2010; Suwa et al. 2010; M¨uller, Janka& Marek 2012; M¨uller, Janka & Heger 2012; Janka et al. 2012; (cid:63) E-mail: [email protected]
Bruenn et al. 2013; Suwa et al. 2013; Pan et al. 2015). However,both the long-time evolution of the successful 2D models as well asthe advent of three-dimensional (3D) simulations have revealed twoserious challenges for the neutrino-driven paradigm: With the ex-ception of Bruenn et al. (2013), the majority of successful 2D sim-ulations tended to produce explosions that are probably underener-getic. Moreover, the most ambitious self-consistent 3D simulationswith multi-group neutrino transport have so far either failed to yieldexplosions (Hanke et al. 2013; Tamborra et al. 2014a,b), or, in thefew successful cases, showed a considerable delay in shock revival(Melson et al. 2015; Lentz et al. 2015) and significantly smaller ex-plosion energies (Takiwaki, Kotake & Suwa 2014). Only the explo-sion of a 9 . M (cid:12) star recently simulated by Melson, Janka & Marek(2015) is an exception from this trend. Whether the neutrino-drivenmechanism provides a robust explanation for shock revival and canaccount for the observed explosion properties of core-collapse su-pernovae may appear doubtful in the light of these results.There is now an emerging consensus about the reasons un-derlying the more fundamental problem of missing or delayed ex-plosions in 3D. Both simple light-bulb and leakage-based simula-tions (Hanke et al. 2012; Couch 2013b,a; Couch & O’Connor 2014;Couch & Ott 2015) as well as models with multi-group transport(Hanke et al. 2013; Takiwaki, Kotake & Suwa 2014) find an ar- c (cid:13) a r X i v : . [ a s t r o - ph . S R ] J u l B. M ¨uller s [ k b / nu c l e on ] s11.0s11.2s11.4s11.60 0.5 1 1.5 2enclosed mass [M O • ] 10 ρ [ g c m - ] s ρ Figure 1.
Density ρ and entropy s for the four progenitors s11.0 (red), s11.2(black), s11.4 (light brown), and s11.6 (green) as a function of enclosedmass. tificial accumulation of turbulent kinetic energy on large spatialscales in 2D due to the action of the inverse turbulent cascade(Kraichnan 1967). E ff ectively, the forward cascade in 3D providesfor more e ffi cient damping / dissipation of the turbulent motions inthe post-shock region, resulting in smaller turbulent kinetic ener-gies. Since the turbulent kinetic energy is directly related to theReynolds stresses (i.e. loosely speaking, the “turbulent pressure”)that have been identified as the primary agent fostering shock ex-pansion in multi-D (Burrows, Hayes & Fryxell 1995; Murphy, Do-lence & Burrows 2013; Couch & Ott 2015; M¨uller & Janka 2015),this a ff ects the critical neutrino luminosity L crit (Burrows & Goshy1993; Murphy & Burrows 2008) required to power an explosiverunaway. However, even though the higher explosion threshold in3D has emerged as a systematic e ff ect, it nonetheless remains asmall e ff ect: Current light-bulb models (Hanke et al. 2012; Couch2013a; Dolence et al. 2013) invariably show a similar reduction of20 . . .
30% in critical luminosity in 2D /
3D compared to 1D, withDolence et al. (2013) still finding a slightly lower explosion thresh-old in 3D. Likewise, neutrino hydrodynamics simulation (Hankeet al. 2013; Takiwaki, Kotake & Suwa 2014) show very similarheating conditions in 2D and 3D prior to shock revival (and eventransient phases with better heating conditions in 3D in Hanke et al.2013), although the small di ff erences between 2D and 3D eventu-ally thwart shock revival in the 3D models of Hanke et al. (2013)and Tamborra et al. (2014a,b). Even though the adverse e ff ects ofthe forward cascade in 3D may still be underestimated by the rela-tively crude grid resolution that current models can a ff ord (Hankeet al. 2012; Couch 2013a; Abdikamalov et al. 2014), it thus emergesthat 3D models of neutrino-driven supernova explosions are veryclose to the explosion threshold. Consequently, relatively small re-finements in the simulations and the initial models may be su ffi -cient to obtain explosions, e.g. moderate rotation (Nakamura et al.2014), magnetic fields (Obergaulinger, Janka & Aloy 2014), as-phericities in the progenitor core (Couch & Ott 2013; Couch et al.2015; M¨uller & Janka 2015), or modifications to the standard neu-trino interaction rates (Melson et al. 2015). The emergence of thespiral mode of the SASI could even allow strongly SASI-dominatedmodels to explode earlier in 3D than in 2D (Fern´andez 2015).By contrast, the problem of underenergetic neutrino-drivenexplosions in 2D has so far gone without a convincing theoret-ical explanation, and fewer suggestions have been made to rem-edy it, although it may be more serious in the sense that it a ff ects even models with successful shock revival: While recent analy-ses of well-observed supernovae (Tanaka et al. 2009; Utrobin &Chugai 2011; Poznanski 2013; Chugai & Utrobin 2014; Pejcha &Prieto 2015) suggest a broader range of explosion energies for typeII-P supernovae between ∼ (0 . . . . × erg instead of a sin-gle “canonical” value of 10 erg, there is arguably still a discrep-ancy, since almost none of the successful 2D and 3D models withmulti-group neutrino transport show explosion energies consider-ably above 10 erg, e.g. the final values at the end of the simula-tions are a few 10 erg in Marek & Janka (2009) for progenitorswith 11 . M (cid:12) and 15 M (cid:12) , (cid:46) erg for a 13 M (cid:12) progenitor in Suwaet al. (2010), and 4 × erg (11 . M (cid:12) progenitor), 1 . × erg(15 M (cid:12) ), and 1 . × erg (27 M (cid:12) ) in Janka et al. (2012). More-over, these estimates are not corrected for the “overburden”, i.e. thebinding energy of the layers outside the shock, so that it remainsunclear whether the explosions become energetic enough to shedthe envelope at all. At first glance, the low explosion energies maysimply be due to the fact that the simulations typically terminatebefore the explosion energy reaches its asymptotic value. While itcan be argued that the final explosion energies may yet be higherby a factor of several because continued accretion sustains strongneutrino emission after shock revival that can power outflows fromthe gain region for (cid:38) . (cid:38) . M grav ≈ . M (cid:12) (Schwab, Podsiadlowski & Rappa-port 2010) inferred from observations (which may, however, su ff erfrom a selection bias). Only the 2D models of Bruenn et al. (2014)form an exception from this trend; they obtain explosion energiesin the range of (3 . . . . . × erg for progenitors between 12 M (cid:12) and 25 M (cid:12) as well as reasonable Nickel masses.While this is encouraging, explosion energies above 10 ergstill remain unexplained, and the problem of underenergetic super-nova explosion will arguably still linger as long as the discrepan-cies between the di ff erent simulation codes are not resolved. In-terestingly, Melson, Janka & Marek (2015) recently reported that3D e ff ects, while apparently detrimental for shock revival in moremassive progenitors, actually boost the explosion energy in a 9 . M (cid:12) progenitor by ∼
10% by reducing the infall velocity in the accretiondownflow and hence the cooling below the gain layer. Here, wefurther investigate their intriguing premise that 3D e ff ects, whilehurtful for shock revival, can prove beneficial in other situationsand contribute to solving the problem of low explosion energies.In this paper, we present a successful 3D multi-group neu-trino hydrodynamics simulation of an 11 . M (cid:12) progenitor with theC o C o N u T-FMT code (M¨uller, Janka & Dimmelmeier 2010; M¨uller& Janka 2015) as further evidence that 3D turbulence plays a posi-tive role after the onset of the explosion. By comparing the dynam-ics and energetics of the 3D explosion model to several long-timesimulations of 2D progenitors in the mass range between 11 M (cid:12) and 11 . M (cid:12) , we demonstrate how 3D e ff ects can lead to a faster,more robust growth of the explosion energy provided that shockrevival can be achieved. So far, the long-time e ff ects of the dimen-sionality during the first hundreds of milliseconds to seconds aftershock revival have received less attention than the role of the thirddimension in shock revival: Successful first-principle models arestill scarce and cannot be evolved for a su ffi ciently long time yetto address this phase in detail, while light-bulb-based studies of su-pernova energetics in 2D and 3D (Handy, Plewa & Odrzywołek2014) cannot adequately account for the feedback of the subsidingaccretion onto the neutrino heating and do not show the drawn-out c (cid:13) , 000–000 ynamics of Neutrino-Driven Supernova Explosions m a x ., a vg ., m i n . s ho c k r a d i u s [ k m ] d i a gno s ti c e xp l o s i on e n e r gy [ e r g ] d E e xp l / d t [ e r g s - ] Figure 2.
Shock propagation and diagnostic explosion energy E expl for the11 . M (cid:12) progenitor in 2D and 3D: The top panel shows the maximum,minimum (solid), and average (dashed) shock radius for model s11.2 2Da(black), s11.2 2Db (blue), and s11.2 3D (red). The middle and bottom panelshow the diagnostic explosion energy E expl and its time derivative d E expl / d t . long-time accretion characteristic of first-principle models. In thispaper, we show that this phase deserves a closer look.Our paper is structured as follows: In Section 2, we reviewthe numerical methods for hydrodynamics and neutrino transportused in our version of the C o C o N u T code, including a brief sketchof the modifications used in the 3D version. In Section 3, we firstgive a descriptive overview of the di ff erences in shock propagationand explosion properties between the 2D and 3D models. Since thequestion of shock revival in 3D is not the objective of our currentstudy, we only provide a brief, cautionary assessment of shock re-vival in the 3D model against the backdrop of recent first-principlemodels in Section 4. In Section 5, we then analyze the physicale ff ects underlying these di ff erences by combining a separate analy-sis of the outflows and downflows in the spirit of Melson, Janka &Marek (2015). Finally, we summarize our results and discuss theirimplications in Section 6. We simulate the collapse and post-bounce phase of the (non-rotating) 11 . M (cid:12) solar-metallicity progenitor s11.2 of Woosley,Heger & Weaver (2002) in 2D and 3D. In order to gauge the ef-fect of stochastic model variations, we perform two di ff erent 2D m a x i m u m s ho c k r a d i u s [ k m ] a v e r a g e s ho c k r a d i u s [ k m ] E e xp l [ e r g ] s11.0_2Ds11.2_2Das11.2_2Dbs11.4_2Ds11.6_2D Figure 3.
Shock propagation and diagnostic explosion energy E expl for thedi ff erent progenitors in 2D: The top and middle panel show the maximumand average shock radius, respectively. The bottom panel shows the diag-nostics explosion energy E expl as a function of time (solid lines). Dashedlines show the time evolution E expl − E ov , i.e. the diagnostic energy cor-rected for the binding energy (overburden) E ov of the material ahead of theshock. Red, black, blue, light brown, and green curves are used for modelss11.0 2D, s11.2 2Da, s11.2 2Db, s11.4 2D, s11.6 2D. simulations (s11.2 2Da and s11.2 2Db ) for this progenitor, andwe also conduct simulations for three other solar metallcity pro-genitors of Woosley, Heger & Weaver (2002) with similar zero-age main sequence (ZAMS) mass and density structure (11 M (cid:12) ,11 . M (cid:12) , 11 . M (cid:12) ). We randomly perturb the radial velocity v r at thebeginning of collapse using a perturbation amplitude δ v r / v r = − .In Figure 1, we show density and entropy profiles for the fourprogenitors simulated in the di ff erent 2D and 3D runs. Despitesmall di ff erences in the location of the interfaces between the dif-ferent shells, the models are very similar in terms of structure with apronounced density jump between the silicon shell and the convec- Scattering on nuclei (including α -particles) was switched o ff after bouncefor model s11.2 2Db, which leads to minor changes in early shock propa-gation, but is inconsequential for the long-time evolution. Since the energyexchange due to recoil in neutrino-nucleon scattering was taken to be pro-portional to the total scattering opacity instead of the neutrino-nucleon scat-tering opacity only (see Equation A31 in M¨uller & Janka 2015) in models11.2 2Da and all other models, the runs with neutrino scattering on nu-clei overestimate pre-heating from heavy flavor neutrinos outside the shockduring the early post-bounce phase (whereas nuclei actually receive negli-gible recoil in scattering reactions), which leads to a slight reduction of theheavy flavor neutrino luminosity and a slightly slower contraction of theproto-neutron star compared to s11.2 2Db.c (cid:13)000
Shock propagation and diagnostic explosion energy E expl for thedi ff erent progenitors in 2D: The top and middle panel show the maximumand average shock radius, respectively. The bottom panel shows the diag-nostics explosion energy E expl as a function of time (solid lines). Dashedlines show the time evolution E expl − E ov , i.e. the diagnostic energy cor-rected for the binding energy (overburden) E ov of the material ahead of theshock. Red, black, blue, light brown, and green curves are used for modelss11.0 2D, s11.2 2Da, s11.2 2Db, s11.4 2D, s11.6 2D. simulations (s11.2 2Da and s11.2 2Db ) for this progenitor, andwe also conduct simulations for three other solar metallcity pro-genitors of Woosley, Heger & Weaver (2002) with similar zero-age main sequence (ZAMS) mass and density structure (11 M (cid:12) ,11 . M (cid:12) , 11 . M (cid:12) ). We randomly perturb the radial velocity v r at thebeginning of collapse using a perturbation amplitude δ v r / v r = − .In Figure 1, we show density and entropy profiles for the fourprogenitors simulated in the di ff erent 2D and 3D runs. Despitesmall di ff erences in the location of the interfaces between the dif-ferent shells, the models are very similar in terms of structure with apronounced density jump between the silicon shell and the convec- Scattering on nuclei (including α -particles) was switched o ff after bouncefor model s11.2 2Db, which leads to minor changes in early shock propa-gation, but is inconsequential for the long-time evolution. Since the energyexchange due to recoil in neutrino-nucleon scattering was taken to be pro-portional to the total scattering opacity instead of the neutrino-nucleon scat-tering opacity only (see Equation A31 in M¨uller & Janka 2015) in models11.2 2Da and all other models, the runs with neutrino scattering on nu-clei overestimate pre-heating from heavy flavor neutrinos outside the shockduring the early post-bounce phase (whereas nuclei actually receive negli-gible recoil in scattering reactions), which leads to a slight reduction of theheavy flavor neutrino luminosity and a slightly slower contraction of theproto-neutron star compared to s11.2 2Db.c (cid:13)000 , 000–000 B. M ¨uller tive shell above the active oxygen burning zone. As we shall see,the 2D models of the di ff erent progenitors are qualitatively verysimilar in terms of explosion dynamics and energetics and shouldthus illustrate the generic behaviour of supernova explosions origi-nating from progenitors in this mass range. The simulations are performed with the general relativistic (GR)neutrino hydrodynamics code C o C o N u T (Dimmelmeier, Font &M¨uller 2002; M¨uller, Janka & Dimmelmeier 2010; M¨uller & Janka2015). Our version of C o C o N u T uses piecewise parabolic (PPM)reconstruction (Colella & Woodward 1984) combined with a hy-brid HLLC / HLLE Riemann solver (Mignone & Bodo 2005) to ob-tain higher-order spatial accuracy. C o C o N u T employs spherical po-lar coordinates ( r , θ, ϕ ), which leads to strong time-step constraintsnear the polar axis in 3D due to the converging grid geometry.We circumvent this problem using an adapted version of the meshcoarsening scheme of Cerd´a-Dur´an (2009): While the equationsof hydrodynamics are solved on a fine grid with constant spac-ing δϕ in longitude everywhere, a filter is applied to the solutionafter each time step to remove short wavelength noise in the ϕ -direction by projecting the conserved variables onto piecewise lin-ear / quadratic functions in “supercells” that contain 2 n ( θ ) fine cellsin the ϕ -direction. The projection algorithm is implemented conser-vatively, and the slopes for the filtered solution are obtained usingthe monotonized-central (MC) limiter of van Leer (1977). The su-percell size 2 n ( θ ) is chosen such that n sin θ > / ∼ ◦ aroundthe pole, which corresponds to 13 .
3% of the total volume. Simi-lar techniques have long been used in numerical meteorology, cp.Kageyama & Sato (2004) and Chapter 18 in Boyd (2001). Polarfiltering allows us to maintain the same e ff ective angular resolutionof 1 . ◦ in 2D and 3D with grids of N r × N θ = ×
128 zones (2D)and N r × N θ × N ϕ = × ×
256 zones (3D, fine grid) coveringthe innermost 10 km of the star, respectively.Like any other solution to avoid the coordinate singularity andthe excessive time step constraint near the axis such as Cartesiancoordinates, overset grids (Kageyama & Sato 2004; Wongwatha-narat, Hammer & M¨uller 2010; Melson, Janka & Marek 2015)or cubed-sphere grids (Ronchi, Iacono & Paolucci 1996; Koldobaet al. 2002; Zink, Schnetter & Tiglio 2008; Fragile et al. 2009),this polar filtering procedure has specific advantages and disadvan-tages: Unlike Cartesian codes, polar filtering allows us to maintainspherical symmetry in the initial conditions and explicit symmetry-breaking terms can be avoided. Di ff erent from cubed-sphere grids,the grid remains orthogonal; and global conservation laws are eas-ier to enforce than for overlapping overset grids. On the other hand,projecting the solution to piecewise linear function e ff ectively in-troduces an anisotropy in the numerical viscosity and di ff usivity(an unwelcome e ff ect that is minimized by overset or cubed-spheregrids but also manifests itself for Cartesian grids that are prone tothe development of m = For example, we use linear functions for the Eulerian density D and themass fractions X i so that the conserved partial masses DX i are representedby quadratic functions. r [ k m ] Figure 7.
Selected mass shell trajectories (black) for model s11.2 3D com-puted from spherically averaged density profiles. The trajectories start withroughly equal spacing in log r shortly before bounce. The plot also showsthe maximum, average, and minimum shock radius (red), the gain radius(light brown), and the radii corresponding to densities of 10 g cm − and10 g cm − (blue). cause the asphericities in the gravitational field are small for non-rotating core-collapse supernovae we use the monopole approxi-mation for the gravitational field, i.e. the lapse function α , the con-formal factor φ , and the radial component β r of the shift vector onlydepend on r , and the non-radial components β θ and β ϕ of the shiftvector are set to zero.For the neutrinos, we use the fast multi-group transport (FMT)scheme of M¨uller & Janka (2015), which is based on a stationarytwo-stream solution of the relativistic transfer equation that is com-bined with an analytic variable Eddington factor closure at low op-tical depths. This schemes includes general relativistic e ff ects un-der the assumption of a stationary metric, but neglects velocity-dependent e ff ects like Doppler shift and aberration. The neutrinorates include emission, absorption, and elastic scattering by nu-clei and free nucleons (along the lines of Bruenn 1985) as well asan e ff ective one-particle rate for nucleon-nucleon bremsstrahlungand an approximate treatment of the energy exchange in neutrino-nucleon scattering reactions. Comparisons of the FMT scheme withthe more sophisticated relativistic neutrino transport solver V er - tex (Rampp & Janka 2002; M¨uller, Janka & Dimmelmeier 2010)showed excellent qualitative and good quantitative agreement. Formore details, we refer the reader to M¨uller & Janka (2015).In order to further alleviate the time-step constraint, the in-nermost part of the computational domain (where densities ex-ceed ∼ × g cm − ) is calculated in spherical symmetry usinga conservative implementation of mixing-length theory for proto-neutron star convection, a procedure that has been used in the con-text of supernova simulations before (e.g. Wilson & Mayle 1988;H¨udepohl 2014). The transition density is adjusted such that it liesinside the convectively stable cooling layer.In the high-density regime, we use the equation of state (EoS)of Lattimer & Swesty (1991) with a bulk incompressibility modu-lus of nuclear matter of K =
220 MeV. At low densities, we employan EoS accounting for photons, electrons and positrons of arbitrarydegeneracy, an ideal gas contribution from baryons (nucleons, pro-tons, α -particles and 14 other nuclear species), Nuclear reactionsare treated using “flashing” as described in Rampp & Janka (2002). c (cid:13) , 000–000 ynamics of Neutrino-Driven Supernova Explosions Figure 4.
Specific entropy for model s11.2 2Da (top row) and model s11.2 3D (in a slice almost perpendicular to the equatorial plane, bottom row) atpost-bounce times of 80 ms, 140 ms, and 181 ms (left to right). Note that a di ff erent color scale for the entropy is used for each of these snapshots. Figure 5.
Specific entropy for model s11.2 2Da (top row) and model s11.2 3D (in a slice almost perpendicular to the equatorial plane, bottom row) atpost-bounce times of 241 ms, 471 ms, and 944 ms (left to right). Note that a di ff erent color scale for the entropy is used for each of these snapshots.c (cid:13)000
Specific entropy for model s11.2 2Da (top row) and model s11.2 3D (in a slice almost perpendicular to the equatorial plane, bottom row) atpost-bounce times of 241 ms, 471 ms, and 944 ms (left to right). Note that a di ff erent color scale for the entropy is used for each of these snapshots.c (cid:13)000 , 000–000 B. M ¨uller
Figure 6.
Volume rendering of the entropy in model s11.2 3D at post-bounce times of 89 ms (top left), 134 ms (top right), 210 ms (bottom left), and 580 ms(bottom right).
Table 1.
Overview of Simulations. The extrapolation of the final remnant masses (last column) is discussed in Section 3.3.Model Progenitor Dimensionality Post-Bounce Diagnostic Energy Baryonic Neutron Star Extrapolated BaryonicTime Reached [s] Reached [erg] Mass Reached [ M (cid:12) ] Remnant Mass [ M (cid:12) ]s11.0 2D s11.0 2 8.195 1 . × . × . × . × . × . × (cid:13) , 000–000 ynamics of Neutrino-Driven Supernova Explosions v s h / v s h , MM s11.0_2Ds11.2_2Das11.2_2Dbs11.4_2Ds11.6_2Ds11.2_3D Figure 8.
Ratio v sh / v sh , MM between the angle-averaged shockvelocity v sh = d r sh , avg / d t and the shock velocity v sh , MM = . E expl / M expl ) / [ M expl / ( ρ pre r )] . predicted by the model ofMatzner & McKee (1999). Note that r sh , avg and its numerical derivativesneed to be smoothed considerably to allow for a useful comparison. v sh / v sh , MM is only shown until 5 s after bounce because the automaticsmoothing procedure becomes ine ff ective toward the end of simulationss11.0 2D, s11.2 2Db, and s11.4 2D and v sh / v sh , MM becomes highlyoscillatory. O • ]-8-6-4-2024 a ng l e - a v e r a g e d r a d i a l v e l o c it y [ c m s - ] Figure 9.
Angle-averaged, density-weighted velocity profiles for models11.2 3D at di ff erent post-bounce times. At the end of the simulation, theangle-averaged velocity is positive outside a mass coordinate of 1 . M (cid:12) ,but the zero point is still moving outward in mass. Note that the angle-average extends over the post-shock and pre-shock region and cannot beused to infer the post-shock velocity directly. In all our simulations, runaway shock expansion sets in when theSi / SiO interface reaches then shock and the mass accretion ratedrops rapidly. Figures 2 (all 2D /
3D 11 . M (cid:12) models) and 3 (long-time evolution of all 2D models) provide an overview over the prop-agation of the shock and the growth of the explosion energies forthe di ff erent models; they show the maximum, minimum (only Fig-ure 2), and angled-averaged shock radius, as well as the “diagnos-tic explosion energy” E expl , which we define as the total net energy(i.e. gravitational + internal + kinetic energy) of all the material thatis nominally unbound and is moving outward with positive radialvelocity at a given time (cp. M¨uller, Janka & Marek 2012; Bruenn M by [ M O • ] s11.0_2Ds11.2_2Das11.2_2Dbs11.4_2Ds11.6_2Ds11.2_3D Figure 10.
Baryonic neutron star masses (comprising all matter at densi-ties higher than 10 g cm − ) for the di ff erent 2D and 3D simulations as afunction of time. et al. 2014). The nucleon rest masses are not included in the internalenergy, i.e. nucleon recombination only contributes to the diagnos-tic energy once it actually takes places. Figure 2 also shows the timederivative of the diagnostic energy. Key results of the simulations,including the diagnostics energy and the baryonic remnant mass atthe end of the simulations as well as estimates for the final remnantmass (see Section 3.3 below), are given in Table 1. ff erences Between 2D and 3D During the First Second For the 11 . M (cid:12) progenitor, the first second after bounce is shownin detail in Figure 2 both in 2D and 3D. In addition, Figures 4and 5 illustrate the multi-dimensional flow morphology for mod-els s11.2 2Da and s11.2 3D on meridional slices, and 3D ray-castimages of neutrino-heated convective bubbles in model s11.2 3Dbefore and after shock revival are shown in Figure 6.Prior to the infall of the Si / SiO interface, we find very similarshock trajectories independent of dimensionality. However, promptconvection develops slightly di ff erently in 2D and 3D, and its resid-ual e ff ect on the entropy and lepton number profiles leads to a slightdivergence between the 2D and 3D models already at early timesin many quantities (neutron star radius, gain radius, cooling pro-files, etc.). This e ff ect is not unphysical per se , but is most probablyexaggerated in our models because the FMT scheme tends to over-estimate the strength of prompt convection. In view of the largesystematic e ff ects that we shall discuss later, it is also inconsequen-tial, but needs to be borne in mind when comparing the di ff erentmodels.After the infall of the Si / SiO interface, the shock expandsslightly faster in 3D than in 2D, and the explosion energy starts toreach appreciable positive values several tens of milliseconds ear-lier. Snapshots of the entropy for models s11.2 3D and s11.2 2Daduring this phase can be seen in the middle and right columns ofFigure 4, which show the development of large convective plumesin both cases. The reader will note that the plumes are somewhataligned with the coordinate axis in 3D, which is clearly a resultof the coordinate choice but need not be considered harmful as dis-cussed in Section 4. At later times, the morphology of the 3D modelis quite di ff erent; instead of the broad, laminar downflows charac-teristic for 2D explosions, the interface between the downflows andthe hot, neutrino-heated ejecta eventually becomes turbulent dur- c (cid:13)000
Baryonic neutron star masses (comprising all matter at densi-ties higher than 10 g cm − ) for the di ff erent 2D and 3D simulations as afunction of time. et al. 2014). The nucleon rest masses are not included in the internalenergy, i.e. nucleon recombination only contributes to the diagnos-tic energy once it actually takes places. Figure 2 also shows the timederivative of the diagnostic energy. Key results of the simulations,including the diagnostics energy and the baryonic remnant mass atthe end of the simulations as well as estimates for the final remnantmass (see Section 3.3 below), are given in Table 1. ff erences Between 2D and 3D During the First Second For the 11 . M (cid:12) progenitor, the first second after bounce is shownin detail in Figure 2 both in 2D and 3D. In addition, Figures 4and 5 illustrate the multi-dimensional flow morphology for mod-els s11.2 2Da and s11.2 3D on meridional slices, and 3D ray-castimages of neutrino-heated convective bubbles in model s11.2 3Dbefore and after shock revival are shown in Figure 6.Prior to the infall of the Si / SiO interface, we find very similarshock trajectories independent of dimensionality. However, promptconvection develops slightly di ff erently in 2D and 3D, and its resid-ual e ff ect on the entropy and lepton number profiles leads to a slightdivergence between the 2D and 3D models already at early timesin many quantities (neutron star radius, gain radius, cooling pro-files, etc.). This e ff ect is not unphysical per se , but is most probablyexaggerated in our models because the FMT scheme tends to over-estimate the strength of prompt convection. In view of the largesystematic e ff ects that we shall discuss later, it is also inconsequen-tial, but needs to be borne in mind when comparing the di ff erentmodels.After the infall of the Si / SiO interface, the shock expandsslightly faster in 3D than in 2D, and the explosion energy starts toreach appreciable positive values several tens of milliseconds ear-lier. Snapshots of the entropy for models s11.2 3D and s11.2 2Daduring this phase can be seen in the middle and right columns ofFigure 4, which show the development of large convective plumesin both cases. The reader will note that the plumes are somewhataligned with the coordinate axis in 3D, which is clearly a resultof the coordinate choice but need not be considered harmful as dis-cussed in Section 4. At later times, the morphology of the 3D modelis quite di ff erent; instead of the broad, laminar downflows charac-teristic for 2D explosions, the interface between the downflows andthe hot, neutrino-heated ejecta eventually becomes turbulent dur- c (cid:13)000 , 000–000 B. M ¨uller ing the infall, resulting in corrugated downflows and partial mixingwith the neutrino-heated ejecta, as can be seen most perspicuouslyin the middle column of Figure 5.Soon after shock revival, the 2D models start to go throughepisodes of halting shock expansion or even transient shock reces-sion. While the growth rate d E expl / d t of the diagnostic explosionenergy reaches values comparable to 3D for 100 . . .
200 ms, the ex-plosion energy grows much less steadily in the long term and hasreached only (4 . . . × erg after 1 s. By contrast, the 3D modelshows a steady growth of the explosion energy (1 . × erg by theend of the simulation), and considerably faster shock expansion. Asillustrated by the mass shell trajectories in Figure 7, the sphericallyaveraged radial velocity behind the shock becomes positives about300 ms after bounce, and the mass shells shocked later than 500 msafter bounce appear to move outward steadily instead of eventuallyfalling back onto the proto-neutron stars. Before we attempt to extrapolate the final remnant masses, it is use-ful to point out a simple analytic relation between the diagnosticenergy and the shock velocity. During the later phases of the explo-sion when the explosion energy has saturated, simple analytic mod-els (Sedov 1959; Kompaneets 1960; Laumbach & Probstein 1969;Klimishin & Gnatyk 1981; Koo & McKee 1990; Matzner & McKee1999) provide a useful qualitative description of shock propagationin hydrodynamical simulations (Woosley & Weaver 1995; Kifoni-dis et al. 2003; Wongwathanarat, M¨uller & Janka 2015). These un-derlying models typically rely on the assumption of self-similarityand / or exponential or power-law approximations for the envelope,neglect the e ff ect of gravity, do not account for continuous energyinput into the ejecta, and have been derived under the assumptionof spherical symmetry. During the first seconds covered in our sim-ulations, all these conditions are violated. It is remarkable that theapproximate formula of Matzner & McKee (1999), v sh , MM = . (cid:114) E expl m (cid:32) m ρ pre r (cid:33) . , (1)nonetheless provides a reasonable estimate for the shock velocity v sh already a few hundreds of milliseconds after shock revival if itis evaluated using appropriate definitions for the explosion energy E expl , the “ejecta” mass m , and the pre-shock density ρ pre : We findthat Equation (1) works well if v sh is taken to be the angle-averagedshock velocity, i.e. the time-derivative of the angle-averaged shockradius r sh , avg , if the pre-shock density is evaluated at r sh , avg , if E expl is identified with the time-dependent diagnostics energy, and if the“ejecta” mass includes the entire mass enclosed by the shock fromabove and the gain radius from below (and thus cannot properly betermed “ejecta” mass as in the original work of Matzner & McKee1999). This is illustrated in Figure 8, which shows the ratio of theangle-averaged shock velocity v sh and the analytic estimate v sh , MM of Matzner & McKee (1999). While there is considerable scatter,the numerical models fall in a band with 1 < v sh / v sh , MM < . v sh = . v sh , MM as a good analytic estimate for earlyshock propagation in core-collapse supernovae. Although the explosion energy in model s11.2 3D has not yetreached its final value and there is still some accretion onto the by +M gain [M O • ]00.511.522.533.54 v s h , MM , β / ( β− ) v e s c β=7 β=4 Figure 11.
Comparison of the post-shock velocity computed from Equa-tion (1) (red curve) to the required shock velocity β/ ( β − v esc for the sep-aration of outgoing and infalling mass shells for two di ff erent values of thecompression ratio β (black curves) for model s11.2 3D. proto-neutron star, there is no doubt that the incipient explosionwill eventually expel the envelope. This is not only clear from thesteady outward movement of the shocked mass shells at late timesvisible in Figure 7 and in the velocity profiles depicted in Figure 9;the diagnostic explosion energy is also significantly higher than theresidual binding energy of the pre-shock matter (the “overburden”in the terminology of Bruenn et al. 2013, 2014) of 5 × erg.This is clearly di ff erent from models s11.2 2Da and s11.2 2Db andthe similar 2D explosion models of the same progenitor discussedin Buras et al. (2006a); Marek & Janka (2009); M¨uller, Janka &Marek (2012).Nonetheless, the long-time simulations of the 2D models overseveral seconds show that even such supposedly tepid models even-tually develop steady shock propagation and reach su ffi ciently highdiagnostic explosion energies to shed the envelope (Figure 1). Forthe low-mass progenitors simulated here, the outgoing and infallingmass shells inevitably separate as the post-shock velocities behindthe entire shock front become positive once it reaches the edge ofthe relatively small C / O core (1 . . . . . M (cid:12) ), which happens al-ready after a few seconds in these models. The acceleration of theshock at the steep density gradient between the C / O core and theHe shell and the small binding energy of the He shell then resultin a steady outward movement of the shocked matter. At that point,the overburden of the unshocked envelope becomes almost negligi-ble (e.g. 10 erg for s11.2 2Db, 5 × erg for s11.6 2D), and wecan determine relatively firm lower limits for the final explosionenergy.Continuous accretion over several seconds provides for suf-ficient neutrino heating to power outflows and pump additionalenergy into the ejecta over this long time-scale, albeit at a rathermodest rate. As a result, models that appear woefully underener-getic during the first second can still develop appreciable explosionenergies, the best example being the 11 . M (cid:12) model, where E expl grows from 3 . × erg at 1 s to 2 . × erg after 11 s. The ex-plosion energies obtained after several seconds are comparable tothe 3D case and compatible with supernova explosion energies atthe lower end of the observed spectrum (see, e.g., Pejcha & Prieto2015).The fact that the diagnostic explosion energies increase moreor less steadily over several seconds in the 2D models (except for c (cid:13) , 000–000 ynamics of Neutrino-Driven Supernova Explosions transient phases where the explosion geometry changes because aneutrino-driven outflow is shut o ff as discussed in Section 5.2.4)has important implications for the usefulness of the diagnostic en-ergy as a predictor of the final explosion properties. Perego et al.(2015) have recently pointed out on the basis of artificial 1D explo-sions of Perego et al. (2015) that the diagnostic energy overshoots the final explosion energy and only approaches its asymptotic valuevery slowly on a time-scale of seconds, and suggest that a better es-timate for the final explosion energy can be obtained by subtractingthe overburden E ov , i.e. the binding energy of the mass shells out-side the shock, from E expl . As illustrated by the comparison of E expl and E expl − E ov in Figure 3, our 2D models show a somewhat dif-ferent behaviour; similar to the simulations of Bruenn et al. (2014),there is no overshooting of E expl above its prospective final valuefor which it appears to furnish a lower bound rather than an upperbound. This is the result of a fundamentally di ff erent way to powerthe explosion in multi-D compared to 1D: Once an explosion istriggered in 1D, the outflow rate quickly drops and only a weakneutrino-driven wind can still pump energy into the ejecta overtime-scales of seconds. The accumulation of shocked material withnegative total energy therefore quickly dominates the total energybudget of the ejecta region and E expl decreases, while E expl − E ov remains roughly constant by virtue of total energy conservation.The case for E expl − E ov as a more compelling predictor of the finalexplosion energy is weaker in multi-D, however, where neutrino-driven outflows can continuously pump energy into the ejecta at ahigh rate, and a considerable part of the shocked material with neg-ative total energy is channeled onto the proto-neutron star insteadof being swept along by the ejecta and reducing the diagnostic en-ergy. E expl may still decrease somewhat on time-scales longer than5 . . .
10 s as the shock propagates through the helium shell, and thisintroduces a residual uncertainty of up to ∼
15% in the final explo-sion energy, which we expect to lie in the range bracketed by E expl and E expl − E ov . It is also noteworthy that E expl − E ov does not ap-pear to be a good predictor for the final explosion energy at earlytimes simply because its rise phase is much more drawn out thanin artificial 1D explosions and it only becomes positive ∼ erg, this comes at the expense of rather high neutronstar masses: Figure 10 shows that the baryonic neutron star masses M by in the 2D models all end up at values (cid:38) . M (cid:12) and will def-initely exceed 1 . M (cid:12) in cases like s11.0 2D and s11.6 2D. Unlessselection e ff ects favor the production of less massive neutron starsin binary systems for some reason, this potentially presents a seri-ous conflict with the inferred neutron star mass distribution. Evenif the lowest-mass neutron stars are presumed to originate fromelectron-capture supernovae, the masses of the neutron stars in the2D models would end up well above the mean value of inferredbaryonic masses of 1 . M (cid:12) (Schwab, Podsiadlowski & Rappaport2010). Since the simulated models represent progenitors with rela-tively small cores and a relatively early onset of the explosion, thisis highly problematic. It is interesting to note that the 2D modelsof Bruenn et al. (2014) also show such a tendency towards highneutron star masses despite their relatively high explosion energies(although this tendency is less striking than in our long-time sim-ulations), with M by = . M (cid:12) for their 12 M (cid:12) model B12-WH07and values well above 1 . M (cid:12) for the three remaining simulations.The faster rise of the explosion energy in 3D could help to re-solve this discrepancy. Although the neutron star mass has not yetconverged to a final value, the spherically-averaged velocity pro-files (Figure 9 indicate that the final “mass cut” is slowly emerging. At the end of the simulation, the net mass accretion rate onto theneutron star in s11.2 3D is lower by a factor of ∼ v post becomes compara-ble to the escape velocity v esc . For a strongly asymmetric explosiongeometry, v post is of course strongly direction-dependent. Hencematerial ahead of the neutrino-heated plumes originating from agiven mass coordinate m in the progenitor will be accelerated to ahigher post-shock velocity by the shock than material with the sameinitial m that is hit later in a direction where the shock expands moreslowly, so that the actual “mass cut” does not correspond to a singlemass shell m in the progenitor. Instead, the dividing line in initialmass coordinate will depend on angle. Nonetheless, one can arguethat the criterion v post = v esc still yields a fairly reliable estimate forthe final mass of the neutron star if an appropriate spherical averagefor v post is used.Equation (1) for the average shock velocity allows us toextrapolate the evolution of the v post if necessary to estimate aspherically-averaged “mass cut”. If the pre-shock velocity is as-sumed to be negligible, the post-shock velocity becomes v post = β − β v sh , (2)in terms of the ratio β of the post- and pre-shock density, and equat-ing this to the escape velocity yields the criterion β − β v sh = (cid:114) G ( M by + M gain ) r , (3)where we include the entire mass interior to the shock and not justthe mass of the neutron star when computing the escape velocity.At late stages, the compression ratio β typically drops below thevalue β = γ = / / or because the strongshock approximation is not strictly applicable over the downflows.We therefore compare v sh in model s11.2 3D to the critical veloc-ity β/ ( β − v esc for two di ff erent values of β in Figure 11. Fig-ure 11 suggests a final baryonic remnant mass of 1 . . . . . M (cid:12) for s11.2 3D (to which late-time fallback might be added). Thiswould imply that the shock has already passed the initial mass cutin some directions. Estimates along the same lines for the long-time simulations in 2D yield baryonic remnant masses of 1 . M (cid:12) for s11.0 2D, 1 . M (cid:12) for s11.4 2D, 1 . M (cid:12) for s11.0 2Db, and1 . M (cid:12) for s11.6 2D assuming β = M grav , M grav ≈ M by − . M (cid:12) (cid:32) M grav M (cid:12) (cid:33) , (4)which provides a reasonable fit across di ff erent nuclear equations ofstates, the estimated baryonic neutron star mass for s11.2 3D can beconverted to a gravitational mass of 1 . M (cid:12) , which would be wellwithin the range of observed neutron star masses and slightly belowthe mean value of the higher-mass population of neutron stars fromiron core progenitors suggested by Schwab, Podsiadlowski & Rap-paport (2010). For the 2D simulations the estimated gravitational Equation (1) is also more convenient to use from the numerical point ofview because the computation of the shock velocity as a numerical deriva-tive of the shock position typically yields very noisy results.c (cid:13)000
15% in the final explo-sion energy, which we expect to lie in the range bracketed by E expl and E expl − E ov . It is also noteworthy that E expl − E ov does not ap-pear to be a good predictor for the final explosion energy at earlytimes simply because its rise phase is much more drawn out thanin artificial 1D explosions and it only becomes positive ∼ erg, this comes at the expense of rather high neutronstar masses: Figure 10 shows that the baryonic neutron star masses M by in the 2D models all end up at values (cid:38) . M (cid:12) and will def-initely exceed 1 . M (cid:12) in cases like s11.0 2D and s11.6 2D. Unlessselection e ff ects favor the production of less massive neutron starsin binary systems for some reason, this potentially presents a seri-ous conflict with the inferred neutron star mass distribution. Evenif the lowest-mass neutron stars are presumed to originate fromelectron-capture supernovae, the masses of the neutron stars in the2D models would end up well above the mean value of inferredbaryonic masses of 1 . M (cid:12) (Schwab, Podsiadlowski & Rappaport2010). Since the simulated models represent progenitors with rela-tively small cores and a relatively early onset of the explosion, thisis highly problematic. It is interesting to note that the 2D modelsof Bruenn et al. (2014) also show such a tendency towards highneutron star masses despite their relatively high explosion energies(although this tendency is less striking than in our long-time sim-ulations), with M by = . M (cid:12) for their 12 M (cid:12) model B12-WH07and values well above 1 . M (cid:12) for the three remaining simulations.The faster rise of the explosion energy in 3D could help to re-solve this discrepancy. Although the neutron star mass has not yetconverged to a final value, the spherically-averaged velocity pro-files (Figure 9 indicate that the final “mass cut” is slowly emerging. At the end of the simulation, the net mass accretion rate onto theneutron star in s11.2 3D is lower by a factor of ∼ v post becomes compara-ble to the escape velocity v esc . For a strongly asymmetric explosiongeometry, v post is of course strongly direction-dependent. Hencematerial ahead of the neutrino-heated plumes originating from agiven mass coordinate m in the progenitor will be accelerated to ahigher post-shock velocity by the shock than material with the sameinitial m that is hit later in a direction where the shock expands moreslowly, so that the actual “mass cut” does not correspond to a singlemass shell m in the progenitor. Instead, the dividing line in initialmass coordinate will depend on angle. Nonetheless, one can arguethat the criterion v post = v esc still yields a fairly reliable estimate forthe final mass of the neutron star if an appropriate spherical averagefor v post is used.Equation (1) for the average shock velocity allows us toextrapolate the evolution of the v post if necessary to estimate aspherically-averaged “mass cut”. If the pre-shock velocity is as-sumed to be negligible, the post-shock velocity becomes v post = β − β v sh , (2)in terms of the ratio β of the post- and pre-shock density, and equat-ing this to the escape velocity yields the criterion β − β v sh = (cid:114) G ( M by + M gain ) r , (3)where we include the entire mass interior to the shock and not justthe mass of the neutron star when computing the escape velocity.At late stages, the compression ratio β typically drops below thevalue β = γ = / / or because the strongshock approximation is not strictly applicable over the downflows.We therefore compare v sh in model s11.2 3D to the critical veloc-ity β/ ( β − v esc for two di ff erent values of β in Figure 11. Fig-ure 11 suggests a final baryonic remnant mass of 1 . . . . . M (cid:12) for s11.2 3D (to which late-time fallback might be added). Thiswould imply that the shock has already passed the initial mass cutin some directions. Estimates along the same lines for the long-time simulations in 2D yield baryonic remnant masses of 1 . M (cid:12) for s11.0 2D, 1 . M (cid:12) for s11.4 2D, 1 . M (cid:12) for s11.0 2Db, and1 . M (cid:12) for s11.6 2D assuming β = M grav , M grav ≈ M by − . M (cid:12) (cid:32) M grav M (cid:12) (cid:33) , (4)which provides a reasonable fit across di ff erent nuclear equations ofstates, the estimated baryonic neutron star mass for s11.2 3D can beconverted to a gravitational mass of 1 . M (cid:12) , which would be wellwithin the range of observed neutron star masses and slightly belowthe mean value of the higher-mass population of neutron stars fromiron core progenitors suggested by Schwab, Podsiadlowski & Rap-paport (2010). For the 2D simulations the estimated gravitational Equation (1) is also more convenient to use from the numerical point ofview because the computation of the shock velocity as a numerical deriva-tive of the shock position typically yields very noisy results.c (cid:13)000 , 000–000 B. M¨uller M ou t · [ M O • s - ] F h , d E e xp l / d t [ e r g s - ] h ou t , e ou t [ M e V / b a r yon ] Figure 12.
Comparison of key properties of the neutrino-driven outflows in2D (black curves) and 3D (red curves) as measured at a radius of 400 km.The top panel shows the mass outflow rate ˙ M out . The middle panel showsthe total enthalpy flux F h , out as defined in Equation (6) (thick lines) along-side the time derivative d E expl / d t of the explosion energy (thin lines). Theaverage total enthalpy ¯ h out (thick lines) and total energy ¯ e out in the outflowsare shown in the bottom panel. masses are higher by more than 0 . M (cid:12) . The 3D e ff ects responsi-ble for the steeper rise of the explosion energy thus improve theagreement with the observational constraints considerably. In the remaining part of the paper, our main thrust will be to ex-plain the physical mechanisms behind the pronounced di ff erencesbetween 2D and 3D simulations in the explosion phase presented inSection 3. We do not investigate the di ff erences in the pre-explosionphase, because the numerical methodology used in this study onlyallows limited conclusions concerning the problem of shock revivalin 3D for reasons detailed below. Nonetheless, a few remarks aboutshock revival in model s11.2 3D are in order, if only to motivatewhy the remainder of this paper focuses completely on the explo-sion phase, and why simulations with a more rigorous treatment ofthe neutrino transport and the neutrino rates are needed to decidethe fate of this particular progenitor model.Superficially, our results for the 11 . M (cid:12) star in 3D may ap-pear to be at odds with the recent core-collapse supernova sim-ulations with multi-group neutrino transport find either no shockrevival at all in 3D or only delayed shock revival compared to 3D(Hanke et al. 2012; Hanke 2014; Tamborra et al. 2014a; Lentz et al. 2015; Melson et al. 2015). Specifically, the 11 . M (cid:12) model failedto explode in 3D (Tamborra et al. 2014a) in a simulation usingthe V ertex -P rometheus code (Rampp & Janka 2002; Buras et al.2006b). However, one should not attach undue importance to thedi ff erent outcomes of the V ertex -P rometheus and C o C o N u T-FMTmodels: Although reasonably close agreement with the more rig-orous transport scheme in V ertex can be reached with the FMTscheme, the di ff erences noted by M¨uller & Janka (2015) are suf-ficiently large to matter for a marginal model like s11.2. The factthat we find an explosion merely underscores how close the 11 . M (cid:12) V ertex model of Hanke (2014) and Tamborra et al. (2014a) comesto an explosive runaway, and that the extreme sensitivity of the su-pernova problem to the neutrino transport treatment and the micro-physics (Lentz et al. 2012b,a; Melson et al. 2015) requires highlyaccurate first-principle models in order to reliably decide whetheran individual progenitor close to the explosion threshold explodesor fails (although “imperfect” simulations may already unearthmuch of the relevant physics from such cases). Considering that thecomparison between the FMT scheme and V ertex revealed slightlybetter heating conditions for a 15 M (cid:12) progenitor at early times, andthat the average shock radius initially expands somewhat further in3D than in 2D in the simulations of the 11 . M (cid:12) progenitor withV ertex , the di ff erent outcome of the FMT and V ertex runs is byno means unexpected.Potentially, the emergence of large-scale bubbles aligned withthe coordinate axis (Figures 4, 5 and 6) also helps in pushing the 3Dmodel above the explosion threshold at an early time (cp. Thomp-son 2000; Dolence et al. 2013; Fern´andez 2015 for the role of thebubble size in the development of runaway shock expansion). Thealignment is clearly a consequence of our coordinate choice andmay also be connected to the coarsening procedure for the polarregion. Such artifacts are unavoidable for standard spherical polar,cylindrical, or Cartesian coordinates (where they manifest them-selves as a preferred excitation of m = ff ective numerical dif-fusivity and viscosity of a code, and physical instabilities will growpreferentially in directions where they are least suppressed (or evenaided) by numerical dissipation. If there is su ffi cient time for insta-bilities like convection and the SASI to reach saturation and gothrough several overturn time-scales or oscillation periods, theseinitial artifacts from the growth phase are eventually washed out,but in a situation where the growth of certain modes acceleratesrapidly (e.g. after the infall of the Si / SiO interface) and then freezesout they can subsist throughout the simulation. Nonetheless, we donot view this a concern; the convective flow does not show any gridalignment prior to the infall of the interface, and outflows eventu-ally develop in the equatorial plane as well. Moreover, we foundno grid alignment of sloshing / spiral motions in SASI-dominatedmodels (which will be reported elsewhere). In more realistic sim-ulations, the explosion geometry will be dictated by anisotropiesin the initial model, e.g. due to convective nuclear burning (Arnett1994; Bazan & Arnett 1994, 1998; Asida & Arnett 2000; Kuhlen,Woosley & Glatzmaier 2003; Meakin & Arnett 2006, 2007b,a; Ar-nett & Meakin 2011; Couch & Ott 2015) or rotation. In a sense,the alignment of the most prominent high-entropy bubbles with theaxis in model s11.2 3D is even fortunate for our further analysisbecause it eliminates the unavoidable grid alignment of 2D explo-sion models as a potential cause for the di ff erent energetics in 2Dand 3D. c (cid:13) , 000–000 ynamics of Neutrino-Driven Supernova Explosions Q · h ea t , Q · c oo l [ e r g s - ] Figure 13.
Volume-integrated heating / cooling rates ˙ Q heat and ˙ Q cool in thegain and cooling region for models s11.2 2Da and s11.2 3D. The innerboundary of the cooling region is defined (somewhat arbitrarily) by a den-sity of 10 g cm − . Note that a di ff erent scale is used for both rates. -70-60-50-40-30-20-100 e t o t (r g a i n ) , − G M / r g a i n [ M e V / b a r yon ] r g a i n [ k m ] T g a i n [ M e V ] Figure 14.
Properties of the gain radius in 2D and 3D: The top panel showsthe binding energy (i.e. the sum of the gravitational, kinetic, and internalenergy) at the gain radius in 2D (black solid curve) and 3D (red solid curve)alongside the Newtonian potential of the neutron star GM / r gain (dashedlines), which sets the typical energy scale at the gain radius. Here, M is theneutron star mass, and r gain is the gain radius. The estimate for e gain fromEquation (13), which is based on the assumption that the Bernoulli integralat the gain radius is zero, is shown in blue for comparison. The middle panelshows r gain itself, and the bottom panel shows the angle-averaged tempera-ture at the gain radius, T gain . η ou t Figure 15.
Outflow e ffi ciency η out , for models s11.2 2Da and s11.2 3D. η out is defined as the ratio between the actual mass outflow rate ˙ M out and afiducial scale ˙ Q heat / | e gain | for the mass loss rate, see Equation (14). Notethat η out is not limited to values η out (cid:54) /
3D DIFFERENCES
We now turn to the underlying physical mechanism responsiblefor the di ff erent evolution of the 2D and 3D models during theexplosion phase. The first step towards understanding the di ff er-ent dynamics of the 2D and 3D models consists in considering theoutflows and downflows separately (in the vein of Melson, Janka& Marek 2015) to analyze the injection of mass and energy intothe “ejecta region” with positive binding energy (similar to Bruennet al. 2014). We partition the computational domain into two re-gions with positive radial velocity ( v r >
0) and negative radial ve-locity ( v r <
0) and then compute total fluxes and averages of severalhydrodynamic quantities. In order not to detract the reader from thephysics, we work with Newtonian definitions here, and the gener-alization to the relativistic case is discussed in Appendix A instead.Unless explicitly stated otherwise, the analysis is based on modelss11.2 2Da and s11.2 3D, i.e. we always refer to model s11.2 2Daand not to model s11.2 2Db when talking about the 2D case.
The first quantities to consider are the mass fluxes ˙ M in and ˙ M out inthe downflows and outflows,˙ M in / out = (cid:90) v r ≶ ρ v r r d Ω , (5)where ρ is the density, and the less and greater signs refer to down-flows and outflows respectively. Furthermore, we define total en-thalpy and energy fluxes F h and F e , F h , in / out = (cid:90) v r ≶ (cid:104) ρ ( (cid:15) + v / + Φ ) + P (cid:105) v r r d Ω (6) The quantity h tot = (cid:15) + P /ρ + v / + Φ , which we shall usually designatein this paper as “total enthalpy” is also referred to as Bernoulli integral orstagnation enthalpy in other contexts (including the case without gravity).In this paper, we prefer the term “total enthalpy” to keep the terminologycompact and stress its close relation to the total energy per unit mass.c (cid:13)000
The first quantities to consider are the mass fluxes ˙ M in and ˙ M out inthe downflows and outflows,˙ M in / out = (cid:90) v r ≶ ρ v r r d Ω , (5)where ρ is the density, and the less and greater signs refer to down-flows and outflows respectively. Furthermore, we define total en-thalpy and energy fluxes F h and F e , F h , in / out = (cid:90) v r ≶ (cid:104) ρ ( (cid:15) + v / + Φ ) + P (cid:105) v r r d Ω (6) The quantity h tot = (cid:15) + P /ρ + v / + Φ , which we shall usually designatein this paper as “total enthalpy” is also referred to as Bernoulli integral orstagnation enthalpy in other contexts (including the case without gravity).In this paper, we prefer the term “total enthalpy” to keep the terminologycompact and stress its close relation to the total energy per unit mass.c (cid:13)000 , 000–000 B. M¨uller F e , in / out = (cid:90) v r ≶ ρ ( (cid:15) + v / + Φ ) v r r d Ω (7)where (cid:15) is the mass-specific internal energy, v is vectorial fluid ve-locity, and Φ is the Newtonian gravitational potential. The rationalefor including the gravitational potential in these fluxes is that thereis a conservation law for the total energy density e tot = (cid:15) + v / +Φ , ∂∂ t (cid:34) ρ (cid:32) (cid:15) + v + Φ (cid:33)(cid:35) + ∇ · (cid:34) ρ (cid:32) (cid:15) + v + Φ (cid:33) v + P v (cid:35) = , (8)if the gravitational potential is time-independent and the conversionof rest mass energy is either disregarded or nuclear rest masses areincluded in the internal energy. Since the diagnostic explosion en-ergy is defined as an integral over the total energy of the ejecta, theenergy budget of the ejecta naturally involves the total enthalpy fluxfrom lower regions of the gain layer to the “ejecta region” with e tot if the boundary of the ejecta region remains at a roughly constantradius.In Figure 12, we show ˙ M out and the average total enthalpyand energy ¯ h tot and ¯ e out (defined as ¯ h out = F h , out / ˙ M out and ¯ e out = F e , out / ˙ M out , and excluding rest mass contributions) in the outflowsat a radius of 400 km as a function of time. This radius has beenchosen because recombination into α -particles, which roughly setsthe final mass-specific total energy in the ejecta (Scheck et al.2006), is already complete at this point, so that F h , out roughlyrepresents the rate at which net total energy is pumped into theejecta region assuming steady-state conditions (i.e. small varia-tions of F h , in with time and radius). This is indeed a very goodapproximation as the comparison for F h , out with the time derivatived E expl / d t of the explosion energy as shown by the middle panelin Figure 12. d E expl / d t correlates extremely well with F h , out , butis slightly smaller. The di ff erence is due to the accumulation ofshocked material with slightly negative total energy and energy ex-change with the downflows due to turbulent di ff usion. The simi-larity of d E expl / d t and F h , out also indicates that explosive nuclearburning does not play a major role for the 11 . M (cid:12) model in agree-ment with earlier 2D simulations with the V ertex -C o C o N u T code(M¨uller, Janka & Marek 2012).On average, the total enthalpy flux into the ejecta region islarger in 3D than in 2D as expected from the di ff erent evolutionof the explosion energy. Interestingly, the relative di ff erence be-tween 2D and 3D in the mass outflow rate ˙ M out is even largerthan for F h , out . The smaller outflow rate in 2D is partially compen-sated by a larger average mass-specific total enthalpy and energyin the outflows, which can be larger than the recombination energy(7 . . . . / nucleon). Thus, care must be exercised in explain-ing di ff erences between 2D and 3D based on the mass outflow rateor the total mass in the gain region alone (Scheck et al. 2006; Mel-son, Janka & Marek 2015) assuming the same contribution to theexplosion energy per unit mass from nucleon recombination into α -particles and heavy nuclei irrespective of the dimensionality. Thatassumption would require similar average enthalpies and energiesin the outflows in 2D and 3D, which is clearly not the case in gen-eral; the di ff erences can be as large as several MeV / nucleon. Re-combination still sets the scale for the asymptotic total energy perunit mass of neutrino-heated ejecta, but hydrodynamic e ff ects mod-ify its precise value in 2D and 3D in di ff erent directions as we shallexplain below.These di ff erences are all the more astonishing because thevolume-integrated neutrino heating rate ˙ Q heat in the gain region(Figure 13) is very similar in 2D and 3D. Especially at late times,˙ Q heat is consistently higher in 2D than in 3D (as is the time-integrated neutrino energy deposition). Assuming that the outflow rate is determined by the total heating rate and the binding energy e gain at the gain radius as ˙ M out ∼ ˙ Q heat / | e gain | , the lower outflow ratein 2D suggests that the material at the gain radius is more stronglybound in this case. This is borne out by Figure 14, which showsthat the binding energy at the gain radius is larger in 2D by up to afactor of ∼ r gain (bottom panel of Figure 14) as a result of whichthe typical energy scale GM / r gain ( M being the neutron star mass)at the gain radius is larger. However, it is evident that this e ff ect can-not fully account for the di ff erence in e gain . The small value of e gain in 3D, indicates that it is not determined by the gravitational energyscale alone. Indeed, a much better estimate for the scale of e gain canbe obtained if we suppose that the Bernoulli integral (including restmasses) at the gain radius is roughly zero, (cid:15) + v + P ρ + Φ = . (9)If we split the internal energy (cid:15) into a thermal component (cid:15) therm and a rest mass contribution (cid:15) rm (normalized to Fe), and assumevanishing velocities as well as an ideal gas equation P = (cid:15) therm / (cid:15) therm + (cid:15) rm + P ρ − GMr gain = , (10)43 (cid:15) therm + (cid:15) rm − GMr gain = , (11) (cid:15) therm = (cid:32) GMr gain − (cid:15) rm (cid:33) , (12)where M is the neutron star mass. For an electron fraction of Y e = . (cid:15) rm would be identical to the recombination energy of (cid:15) rec ≈ . / nucleon from protons and neutrons with equal massfractions into iron group elements, and for our purposes this stillprovides a su ffi cient approximation even though we have Y e < . e gain (in which we ex-cluded rest masses), we thus obtain e gain ≈ (cid:15) therm − GMr gain ≈ − (cid:15) rec − GM r gain . (13)As shown in Figure 14, this still overestimates | e gain | a bit, but ac-counts for the slow rise of | e gain | compared to the gravitational en-ergy scale GM / r gain during the contraction of the neutron star.While the di ff erent absolute value of the binding energy atthe gain radius is part of the explanation for the di ff erent massoutflow rates, there is also an additional e ff ect at play. In gen-eral, the mass outflow rate will only be approximately given by˙ M out ∼ ˙ Q heat / | e gain | , and one can introduce an e ffi ciency parameter η out to compare the actual mass outflow rate with this fiducial rate, η out = | e gain | ˙ M out ˙ Q heat . (14) η out is plotted in Figure 15. On average, the outflow e ffi ciency η out isalso considerably larger in 3D (where it fluctuates around η out ≈ η out ≈ . e gain at the gain radius, the conversion of neutrinoheating into an outflow is less e ffi cient. c (cid:13) , 000–000 ynamics of Neutrino-Driven Supernova Explosions Figure 16.
Snapshots of the specific entropy s (left column, measured in k b / nucleon) and the radial velocity v r (right column, measured in cm s − ) for models11.2 2Da (left halves of the individual panels) and for a slice of model s11.2 3D (right halves) at a post-bounce time of 400 ms. The same data is shown in allplots, only the zoom level is di ff erent. Note the broad equatorial accretion downflow and the formation of a secondary accretion shock at a radius of ∼
100 kmin 2D.
These two e ff ects, as well as the higher asymptotic energy per unitmass in 2D can be traced to the constrained axisymmetric flow ge-ometry and a fundamentally di ff erent behaviour of the accretiondownflows in 2D compared to 3D. The di ff erent flow morphologyis illustrated qualitatively in Figure 16, which shows snapshots (for di ff erent spatial scales) of the radial velocity and entropy in 2D and3D for a representative post-bounce time of 400 ms. c (cid:13)000
These two e ff ects, as well as the higher asymptotic energy per unitmass in 2D can be traced to the constrained axisymmetric flow ge-ometry and a fundamentally di ff erent behaviour of the accretiondownflows in 2D compared to 3D. The di ff erent flow morphologyis illustrated qualitatively in Figure 16, which shows snapshots (for di ff erent spatial scales) of the radial velocity and entropy in 2D and3D for a representative post-bounce time of 400 ms. c (cid:13)000 , 000–000 B. M¨uller
These snapshots reveal that the interface between the outflowsand the colder, low-entropy becomes turbulent due to the Kelvin-Helmholtz instability in 3D, which distorts the downflows as theyapproach the proto-neutron star (as already mentioned in Sec-tion 3), whereas Kelvin-Helmholtz instabilities are noticeably ab-sent at the shear interfaces between the broad equatorial downflowand the polar bubbles in 2D (middle and bottom row in Figure 16).The tendency of the downflows to become more turbulent in 3Dhas been recognized already by Melson, Janka & Marek (2015) intheir 9 . M (cid:12) model, although the morphological di ff erence between2D and 3D is much more pronounced in a continuously accretingmodel like the 11 . M (cid:12) progenitor. Moreover, although there maybe a deeper connection between the two phenomena, the stabilityof the shear interfaces in 2D is likely due to a di ff erent reason thanthe emergence of large-scale structures in the pre-explosion phase(Hanke et al. 2012) that has been explained by the inverse turbulentenergy cascade in 2D (Kraichnan 1967). Despite the inverse turbu-lent cascade, subsonic shear layers / interfaces remain prone to theKelvin-Helmholtz instability in 2D; and 2D supernova simulationsare easily able to resolve the instability (M¨uller, Janka & Heger2012; Fern´andez 2015) even without extraordinary high resolution.The situation changes, however, for the supersonic shear in-terfaces between the downflows and the neutrino-heated bubblesthat we encounter during the explosion phase. Here, the classi-cal growth rate ω = k ∆ u / k is the wave vector, and ∆ u is the transverse velocity jump across the interface) in the vortexsheet approximation for incompressible flow is no longer applica-ble. Instead modes with a su ffi ciently high e ff ective Mach numberMa = ∆ u cos θ/ c s (where c s is the sound speed and θ is the anglebetween the wave vector and the vectorial velocity jump) are stabi-lized (Gerwin 1968), although the stability analysis is more com-plicated if finite-width shear layers are considered (Blumen 1970;Blumen, Drazin & Billings 1975; Drazin & Davey 1977; Choud-hury & Lovelace 1984; Balsa & Goldstein 1990). This implies thatthe Kelvin-Helmholtz instability can be partially or completely sup-pressed in 2D (where cos θ = θ can be arbitrarily small. In principle, it isconceivable that numerical di ff usivity and viscosity further help tosuppress the Kelvin-Helmholtz instability more strongly and ear-lier than the physics might dictate, but the fact that the instabil-ity evidently operates in the 3D model provides evidence that thenumerical resolution cannot be faulted for the behavior of the 2D It is noteworthy that laser-driven plasma experiments may be able tocapture this e ff ect and quantify the reduced growth or suppression of theKelvin-Helmholtz instability in 2D (Malamud et al. 2013). Loosely speaking, a small value of cos θ guarantees that sound waveson either side of the vortex sheet can propagate in both directions alongthe wave vector k of a given perturbation mode to mediate the pressurefeedback required for the growth of the Kelvin-Helmholtz instability: Fora given vectorial velocity ± u / n and velocity c s in either of the fluids willhave a velocity component ( ± u / + n c s ) · k / | k | = ± u / θ + c s n · k / | k | along the direction of k in the rest frame. If cos θ is su ffi ciently small, thisvelocity component can take on either sign depending on n in both fluids. In2D, we always have ± u / · k / | k | = u /
2, and sound waves cannot propagatein both directions any longer for su ffi ciently large u . Note that this is only aheuristic explanation that cannot predict the critical Mach number correctly;Section B in Gerwin (1968) and the other aforementioned references shouldbe consulted for a rigorous derivation of the dispersion relation. models. Furthermore, the stability of the accretion downflows is apersistent feature even in high-resolution 2D models with continu-ing accretion (Bruenn et al. 2014); it is thus without doubt physicalin origin.The “turbulent braking” of the downflows in 3D is reflectedquantitatively in radial profiles of the average velocity ¯ v in / out , en-tropy ¯ s in / out and mass-specific total energy ¯ e tot , rm , in / out of the down-flows and outflows. We define these quantities as density-weightedaverages (denoted by bars) of the radial velocity v r , the specificentropy s , and the total energy e tot as¯ v in / out = (cid:82) v r ≶ ρ v r d Ω (cid:82) v r ≶ ρ d Ω (15)¯ s in / out = (cid:82) v r ≶ ρ s d Ω (cid:82) v r ≶ ρ d Ω (16)¯ e tot , in / out = (cid:82) v r ≶ ρ e tot , rm d Ω (cid:82) v r ≶ ρ d Ω , (17)and show the results for models s11.2 2Da and s11.2 3D at a post-bounce time of 400 ms in Figure 17. Moreover, we consider radialprofiles of the mass and total enthalpy fluxes ˙ M in / out and ˙ F h , rm , in / out in both streams in Figure 18; these are computed according toEquation (5) and (6). However, for computing radial profiles weinclude rest mass contributions in the total energy and the total en-thalpy flux (as denoted by the additional subscript rm). This defini-tion has the advantage that both ˙ M in , out and ˙ F h , rm , in / out are constantin the limit of stationary streams without mass, energy, and mo-mentum exchange, so that changes in these fluxes serve as usefulindicators for lateral mixing between the outflows and downflows.Due to turbulent braking (i.e. by an e ff ective turbulent eddyviscosity), the average infall velocity in the downflows reaches only1 . × cm s − in 3D, and decreases in magnitude once the down-flows penetrate further down than a radius of ∼
200 km. By contrast,the downflows reach a sizable fraction of the free-fall velocity in 2Dbefore they are abruptly decelerated at a secondary accretion shockat r ≈
100 km. However, the outflow velocities are also higher in2D. During the phase considered here, the entropy of the down-flows does not vary considerably in 2D between r ≈
100 km and r ≈ F h , rm , in / out / ˙ M in / out (bottom panel of Figure 18) does not change ap-preciably in this region. This is a further indication for a lack of dis-sipation by turbulent eddy viscosity and of lateral mixing betweenthe downflows and outflows. By contrast, the slope in ¯ s in and ¯ s out and the “flux-averaged” total enthalpy F h , rm , in / out / ˙ M in / out points tolateral mixing between the downflows and outflows in 3D. The in-crease of F h , rm , in / out / ˙ M in / out in the downflows during the infall fromthe shock is mirrored by a decrease of F h , rm , in / out / ˙ M in / out with radiusin the ouflows due to turbulent mixing of the hot ejecta with coldermatter from the downflows: This explains why the neutrino-heatedejecta contribute only ∼ . . . / nucleon to the diagnostic ex-plosion energy, instead of the 7 . . . . / nucleon available fromnucleon recombination.There is thus ample evidence that turbulent viscosity and dif-fusion brake the accretion funnels and transfer energy from the out-flows. It is tempting to invoke this as an explanation for the lowerbinding energy at the gain radius in 3D (Figure 14). At first glance, c (cid:13) , 000–000 ynamics of Neutrino-Driven Supernova Explosions the fact that both the downflows and outflows are less stronglybound in 3D at the bottom of the gain layer (bottom panel of Fig-ure 17) may seem to conflict with this assumptions, but the lowerbinding energies of the outflows at small radii is to be expected be-cause the supply for outflow comes from freshly accreted matterthat has undergone turbulent braking in the downflows. Turbulentbraking and turbulent di ff usion are perfectly acceptable explana-tions for internal energy distribution within the gain region (butnot the higher enthalpy flux into the ejecta region, see below), andmay thus account for the lower | e gain | in 3D and hence for the highermass outflow rates.Moreover, the turbulent braking in 3D may have implica-tions for the final neutron star masses. In Section 3.3, we assumedthat the “mass cut” occurs roughly when the shock accelerates thenewly swept-up material to the escape velocity. Without e ffi cientlateral momentum transfer and without pressure support from anexpanding hot bubble from below, it seems inevitable that materialover existing downflows must eventually fall onto the neutron starif this condition is not met. Since the free-fall time scale at radiiof several thousands of kilometers (where the initial mass cut esti-mated in Section 3.3 is located) is of the order of seconds, accretionmust necessarily last over a correspondingly long duration. On theother hand, if the downflows are braked by a turbulent eddy vis-cosity below a certain radius in 3D, slowly moving matter in thedownflows may instead be entrained by the high-entropy bubblesin regions where the angular-averaged velocity is already positive,so that the residual accretion onto the neutron star may cease ear-lier and the total mass accreted during the explosion phase may beconsiderably lower than estimated in Section 3.3. Longer 3D sim-ulations will be necessary to investigate this hypothesis.However, the internal redistribution of energy within the gainregion in 3D cannot account for the significantly higher total en-thalpy flux in the outflows and the faster rise of the explosion en-ergy: The higher outflow rate will come at the expense of energyloss from the outflows to the downflows – the overall conserva-tion law cannot be cheated. Consequently, there must be additionalmechanisms that remove energy from the gain region in 2D and re-duce the outflow e ffi ciency η out (Figure 15) compared to 3D. Thedi ff erent dynamics of the outflows and downflows nonetheless re-mains a crucial element of the explanation because it provides thebasis for three mechanisms discussed in the subsequent sections. The lack of turbulent braking in 2D implies that the accretion fun-nels either hit the neutron star directly with a high impact velocityor are decelerated abruptly in a secondary accretion shock (top rowof Figure 16). Even in the latter case, thin accretion funnels stillpenetrate the hot, neutrino-heated matter all the way down to thegain radius, cutting o ff a confined high-entropy bubble from theoutflows. In the snapshots shown in Figures 16 and 19 (with aneven higher zoom level), these narrow downflows strike the proto-neutron star surface with velocities of up to 6 × cm s − . As aresult they overshoot considerably into the convectively stable cool-ing layer and excite strong wave activity. The emission of strongacoustic waves that steepen into shocks and then dissipate is imme-diately evident from the top right panel of Figure 16, but the decel-eration of the downflows also excites strong g -modes in the neu-tron star surface region (a phenomenon that has been thoroughlyanalyzed in the context of gravitational wave emission, cp. Marek,Janka & M¨uller 2009; Murphy, Ott & Burrows 2009; M¨uller, Janka& Marek 2013). In 3D, overshooting is much less pronounced, and -4-3-2-1012 v i n , v ou t [ c m s - ] s i n , s ou t [ k b / nu c l e on ]
30 100 1000 4000r [km]-15-10-505 e t o t , i n , e t o t , ou t [ M e V / nu c l e on ] Figure 17.
Radial profiles of the average velocity (top panel), entropy (mid-dle panel), and total energy (bottom) per nucleon in the outflows (thin lines)and downflows (thick lines) in 2D (black) and 3D (red) at a post-bouncetime of 400 ms. Note that rest mass contributions are included in the totalenergy here. so is the excitation of acoustic waves and g -modes. This can be seenfrom the conspicuous absence of sawtooth-like features in the ve-locity and by considering the radial velocity dispersion (cid:104) ( v r −(cid:104) v r (cid:105) ) (cid:105) ,which is significantly smaller in 3D below the gain radius (toppanel of Figure 20). This result is in agreement with other 3Dsimulations of supernova explosions using self-consistent neutrinotransport (Melson, Janka & Marek 2015) and parameterized neu-trino heating (Murphy, Dolence & Burrows 2013; Handy, Plewa &Odrzywołek 2014) and linear theory, which suggests that the exci-tation of waves ( g -modes in particular) at convective boundaries isstrongly sensitive to the convective Mach number Ma and becomesvery e ffi cient for Ma ∼ g -modes constitutes a non-advective energy drain in 2D; it trans-ports energy from the lower layers of the gain regions deep into thecooling region without the need to transport mass. If the g -modeenergy flux is su ffi ciently high, it provides a very likely explana-tion for the permanently higher binding energy at the gain radius in2D. Unfortunately, the g -mode energy flux in our simulations can-not readily be quantified; this would not only require performinga full spherical Reynolds decomposition, but also detailed knowl-edge about the dispersion relation of the g -modes, which is beyondthe scope of this paper. However, since the transfer of the kineticenergy from the downflows into g -modes involves P d V -work ontothe neutron star surface and any turbulent energy flux into deeperlayers should show up in correlated pressure and velocity fluctua- c (cid:13)000
Radial profiles of the average velocity (top panel), entropy (mid-dle panel), and total energy (bottom) per nucleon in the outflows (thin lines)and downflows (thick lines) in 2D (black) and 3D (red) at a post-bouncetime of 400 ms. Note that rest mass contributions are included in the totalenergy here. so is the excitation of acoustic waves and g -modes. This can be seenfrom the conspicuous absence of sawtooth-like features in the ve-locity and by considering the radial velocity dispersion (cid:104) ( v r −(cid:104) v r (cid:105) ) (cid:105) ,which is significantly smaller in 3D below the gain radius (toppanel of Figure 20). This result is in agreement with other 3Dsimulations of supernova explosions using self-consistent neutrinotransport (Melson, Janka & Marek 2015) and parameterized neu-trino heating (Murphy, Dolence & Burrows 2013; Handy, Plewa &Odrzywołek 2014) and linear theory, which suggests that the exci-tation of waves ( g -modes in particular) at convective boundaries isstrongly sensitive to the convective Mach number Ma and becomesvery e ffi cient for Ma ∼ g -modes constitutes a non-advective energy drain in 2D; it trans-ports energy from the lower layers of the gain regions deep into thecooling region without the need to transport mass. If the g -modeenergy flux is su ffi ciently high, it provides a very likely explana-tion for the permanently higher binding energy at the gain radius in2D. Unfortunately, the g -mode energy flux in our simulations can-not readily be quantified; this would not only require performinga full spherical Reynolds decomposition, but also detailed knowl-edge about the dispersion relation of the g -modes, which is beyondthe scope of this paper. However, since the transfer of the kineticenergy from the downflows into g -modes involves P d V -work ontothe neutron star surface and any turbulent energy flux into deeperlayers should show up in correlated pressure and velocity fluctua- c (cid:13)000 , 000–000 B. M¨uller -6-4-20246 F h , i m , n / ou t [ e r g s - ] -0.1-0.0500.050.1 M i n / ou t · [ M O • s - ]
30 100 1000 4000r [km]-20246810 F h , r m , i n / ou t / M ou t · [ M e V / nu c l e on ] Figure 18.
Overview of the mass and energy fluxes in the outflows anddownflows in 2D. The top panel shows the total enthalpy fluxes F h , rm , in / out in the outflows (positive values) and downflows (negative values) in 2D(black) and 3D (red) at a post-bounce time of 400 ms. The middle panelshows the mass inflow / outflow rates ˙ M in / out , and the bottom panel showsthe average flux-weighted total enthalpies F h , rm , in / out / ˙ M in / out . Note that restmass contributions are included in the total entropy here, unlike in Fig-ures 12 and 17. tions in a transition layer between the convective gain region andthe convectively stabilized cooling layer, we can formulate a crudeestimate for the g -mode flux by computing what is nominally anacoustic luminosity, namely, L P d V = r (cid:90) δ P δ v r d Ω , (18)where δ P and δ v r denote the deviations of the pressure and radialvelocity from their respective angular averages. The resulting esti-mates for the flux are shown in the bottom panel of Figure 20 andpoint to a sizable energy flux of the order of several 10 erg s − from the vicinity of the gain radius into the deeper layers of theproto-neutron star surface (carried by g -modes) and to the outerregions of the gain layer (carried by acoustic waves). Such largefluxes are comparable to the typical total enthalpy flux in the out-flows and even to the volume-integrated neutrino heating rate, andtherefore need to be accounted for in the total energy budget of thegain region.Incidentally, the excitation of acoustic waves also providesan explanation for the high entropy (Figure 17) and total enthalpy(Figures 12 and 18) in the outflows in 2D, which can still increasesomewhat at radii where neutrino heating is essentially irrelevant.The dissipation of the the acoustic waves in the expanding hot bub-ble helps to increase the total energy and entropy in the ejecta re-gion beyond the ∼ . . . . ff ect cannot compensate for the lower mass outflowrate in 2D, and the net e ff ect of wave excitation in 2D remains adetrimental one.The importance of wave excitation at the convective bound-ary will inevitably vary between di ff erent 2D models depending onthe explosion geometry and the duration of accretion. If neutrinoheating is strong, the explosion energy rises steeply, and shock ex-pansion is fast as in the models of Bruenn et al. (2013, 2014), itmay play a less prominent role. The formation of secondary accre-tion shocks as in our 2D models is likely to increase the energy lossby wave excitation tremendously because the e ffi ciency of this pro-cess also depends on the frequency overlap between the convectiveforcing and the excited modes (Goldreich & Kumar 1990; Lecoanet& Quataert 2013). The formation of a secondary accretion shockprovides for fluctuations with typical frequencies inversely propor-tional to the short sound-crossing time-scale (of the order of mil-liseconds) in the confined bubble. Furthermore, the stochastic forc-ing of g -modes in 2D by one or two strong downflows is presum-ably also more e ffi cient than in 3D, where there are more smallerand uncorrelated downflows.Finally, we comment on similarities and di ff erences between g -mode and acoustic wave excitation between our models andthe acoustically-driven explosion models of Burrows et al. (2006,2007), where a strong flux acoustic waves, excited by an (cid:96) = core g -mode with amplitudes of several kilometers, is responsiblefor shock expansion in the first place. Our 2D models are similar tothose of Burrows et al. (2006, 2007) only in the sense that energydeposition by acoustic waves contributes to the growth of the explo-sion energy, but di ff erent from their simulations the acoustic con-tributions remains subdominant compared to the volume-integratedneutrino heating rate, which is more than four times larger at thetime shown in Figure 20 (a constellation that Burrows et al. 2007anticipated when postulating a “hybrid mechanism” with combinedheating by neutrinos and acoustic waves). Moreover, the excitationmechanism for acoustic waves is genuinely di ff erent in our case;they are excited directly by the interaction of the downflows withthe convectively stable neutron star surface layer without the needto channel the accretion power through an (cid:96) = g -mode as a“transducer” as in the models of Burrows et al. (2006, 2007). Sucha large-amplitude core g -mode is not found in our simulations (andcould not have arisen simply because of the spherically symmet-ric treatment of the neutron star core), and the outer g -modes ofrather modest amplitude excited in our models could not act as ane ffi cient transducer in the vein of Burrows et al. (2006, 2007) dueto neutrino losses (see below). It is interesting to note that directexcitation at the convective boundary still allows acoustic waves tocontribute (albeit at a minor level) to the explosion energy in 2Deven without such a transducer.Acoustic energy deposition also remains a secondary e ff ectinsofar as this direct excitation mechanism works e ffi ciently only after the onset of the explosion once the typical Mach number ofthe downflows is su ffi ciently high. Moreover, contrary to Burrowset al. (2006, 2007) the net e ff ect of wave excitation in our models isstill harmful because the power pumped into outer g -modes consti-tutes an energy drain that outweighs the rate of energy depositionby acoustic waves by far. Di ff erent from their models where theenergy in the core g -mode is eventually “recycled” into an acousticenergy flux that drives shock expansion, the energy pumped into theouter g -mode is manifestly lost due to neutrino cooling in our case,and the mode coupling analysis of Weinberg & Quataert (2008) c (cid:13) , 000–000 ynamics of Neutrino-Driven Supernova Explosions suggests that this should also happen if the core g -mode exciteddue to non-linear mode coupling. In addition to energy loss by wave excitation, which contributes tothe higher binding energy at the gain radius, the growth of the ex-plosion energy in 2D is further hampered by the fact that much ofthe neutrino energy deposition occurs in regions where the heatedmatter cannot directly escape in an outflow, i.e. either directly inthe accretion funnels or in high-entropy bubbles confined by down-flows and a secondary accretion shock like the equatorial bubble inFigures 16 and 19, a phenomenon for which we borrow the term“steric hindrance” from chemistry. In principle, such bubbles couldeventually push the secondary accretion shock out by undergoing“secondary shock revival”, but as long as the amount of heating isinsu ffi cient, the bubble cannot break through the surrounding andoverlying downflows.In the snapshot shown in the right panel of 19, the surfacefraction covered by the confined bubble and the downflows exceeds50%, and the heating rate per unit mass is also largest in the down-flows. Since the fast downflows generally occupy a surface fractionof (cid:38)
50% in 2D outside the typical location of a secondary shockat (cid:46) ffi ciency oscillates around η out ∼ . ffi ciencyis of order η out ∼ Finally, the outflows in axisymmetric 2D simulations are less “sta-ble” than in 3D in other respects as illustrated in Figure 22: Whilethe Kelvin-Helmholtz instability between the downflows and out-flows is largely suppressed in 2D, this also implies that plumesof cold material can penetrate far into the neutrino-heated high-entropy bubbles provided that they develop in the first place. Be-cause of the symmetry of the system, these plumes are actuallytoroidal structures, and can therefore completely constrict an out-flow if they reach the symmetry axis (Figure 22). Similarly, a down-flow that wanders toward the pole can also constrict a polar outflowand cut it o ff from fresh supply of neutrino-heated material.These events typically reduce the surface fraction covered byoutflows at the recombination radius (where they start to contributeto the diagnostics explosion energy) for a considerable amount oftime and thus reduce the rate of increase of E expl . Very often theexplosion geometry changes dramatically after such an event andthe surface fraction of the outflows remains small permanently. Insome cases, a high-entropy bubbles is cut o ff completely from thesupply of neutrino-heated matter from below (Figure 23) for sev-eral seconds. Even if an outflow is eventually reestablished in thesame direction, or if it is strong enough to survive because the coldplumes reach the axis at a relatively large radius (Figure 22), theejecta will then typically contain a large amount of cold materialwhose total net energy is barely positive, and the growth of the ex-plosion energy will still be delayed. Such events explain excursionsor even a permanent drop of the the average total enthalpy ¯ h tot inthe outflows to low values < . ff erent guise and occursonly episodically, but with a more catastrophic e ff ect than in 3D.Interestingly, there even appears to be an e ff ect that compensatessomewhat for the suppression of the Kelvin-Helmholtz instabilitiesin 2D due to the supersonic velocities in the downflows: Rayleigh-Taylor instabilities between the high-entropy bubbles and the coldoverlying post-shock matter develop more readily in 2D. This is anatural consequence of higher entropies in the neutrino-heated bub-bles in 2D (middle panel of Figure 17), which implies a higher At-wood number between the bubbles and the colder post-shock mat-ter. Thus, the lack of mixing by Kelvin-Helmholtz instabilities in2D and the entropy boost due the dissipation of acoustic waves canalso have a detrimental side e ff ect on the robustness of the explo-sion. In the most extreme cases of outflow constriction in 2D, the out-flows are shut o ff altogether, and the outflow surface fraction dropsto zero permanently, or at least over several seconds (bottom panelof Figure 21 and Figure 23). This does not imply, however, thatthe explosion has failed; it only implies that neutrino heating isnot strong enough to establish a wind that prevents the fallback ofslowly moving matter in the wake of the shock. The pockets ofcold, slowly moving matter from the C / O shell in the 2D modelsthat will undergo this kind of “early fallback” (bottom right panelof Figure 23) only contain a few hundredths of a solar mass by theend of the simulations, and therefore will not change the neutronstar mass considerably. Moreover, the mass accretion rate onto thesecondary accretion shock is so low at late times that it can startto expand after “secondary shock revival”, thus re-establishing anoutflow (bottom right panel of Figure 23).While not indicative of a failure of the explosion, the smallor vanishing outflow surface fraction in the long-time simulationsnonetheless indicates (like the models of Bruenn et al. 2014) thatthe separation of outgoing and infalling mass shells in 2D worksdi ff erently from the usual picture where a high-entropy neutrino-driven wind with an approximately spherical flow geometry even-tually develops. The polar outflows can be viewed as a confinedwind driven jointly by neutrino heating and acoustic waves, butthey never cover the entire sphere, and because of their flow geom-etry and the strong activity of acoustic waves, the outflow dynamicsis completely di ff erent from spherical winds driven purely by neu-trino heating.From the foregoing, it is clear that the inhibition of theneutrino-driven wind in 2D is probably largely artificial; and weonly mention this peculiarity for that reason. As discussed in Sec-tion 5.2.1, the presence of a larger e ff ective eddy viscosity in 3Dcould terminate accretion earlier than in 2D (where supersonicallyinfalling matter is hardly decelerated by lateral momentum trans-fer), or at least decelarate infalling matter su ffi ciently to be sweptalong by an incipient spherical wind after a few seconds. Moreover,our models likely underestimate the di ff usive neutrino luminosity c (cid:13)000
50% in 2D outside the typical location of a secondary shockat (cid:46) ffi ciency oscillates around η out ∼ . ffi ciencyis of order η out ∼ Finally, the outflows in axisymmetric 2D simulations are less “sta-ble” than in 3D in other respects as illustrated in Figure 22: Whilethe Kelvin-Helmholtz instability between the downflows and out-flows is largely suppressed in 2D, this also implies that plumesof cold material can penetrate far into the neutrino-heated high-entropy bubbles provided that they develop in the first place. Be-cause of the symmetry of the system, these plumes are actuallytoroidal structures, and can therefore completely constrict an out-flow if they reach the symmetry axis (Figure 22). Similarly, a down-flow that wanders toward the pole can also constrict a polar outflowand cut it o ff from fresh supply of neutrino-heated material.These events typically reduce the surface fraction covered byoutflows at the recombination radius (where they start to contributeto the diagnostics explosion energy) for a considerable amount oftime and thus reduce the rate of increase of E expl . Very often theexplosion geometry changes dramatically after such an event andthe surface fraction of the outflows remains small permanently. Insome cases, a high-entropy bubbles is cut o ff completely from thesupply of neutrino-heated matter from below (Figure 23) for sev-eral seconds. Even if an outflow is eventually reestablished in thesame direction, or if it is strong enough to survive because the coldplumes reach the axis at a relatively large radius (Figure 22), theejecta will then typically contain a large amount of cold materialwhose total net energy is barely positive, and the growth of the ex-plosion energy will still be delayed. Such events explain excursionsor even a permanent drop of the the average total enthalpy ¯ h tot inthe outflows to low values < . ff erent guise and occursonly episodically, but with a more catastrophic e ff ect than in 3D.Interestingly, there even appears to be an e ff ect that compensatessomewhat for the suppression of the Kelvin-Helmholtz instabilitiesin 2D due to the supersonic velocities in the downflows: Rayleigh-Taylor instabilities between the high-entropy bubbles and the coldoverlying post-shock matter develop more readily in 2D. This is anatural consequence of higher entropies in the neutrino-heated bub-bles in 2D (middle panel of Figure 17), which implies a higher At-wood number between the bubbles and the colder post-shock mat-ter. Thus, the lack of mixing by Kelvin-Helmholtz instabilities in2D and the entropy boost due the dissipation of acoustic waves canalso have a detrimental side e ff ect on the robustness of the explo-sion. In the most extreme cases of outflow constriction in 2D, the out-flows are shut o ff altogether, and the outflow surface fraction dropsto zero permanently, or at least over several seconds (bottom panelof Figure 21 and Figure 23). This does not imply, however, thatthe explosion has failed; it only implies that neutrino heating isnot strong enough to establish a wind that prevents the fallback ofslowly moving matter in the wake of the shock. The pockets ofcold, slowly moving matter from the C / O shell in the 2D modelsthat will undergo this kind of “early fallback” (bottom right panelof Figure 23) only contain a few hundredths of a solar mass by theend of the simulations, and therefore will not change the neutronstar mass considerably. Moreover, the mass accretion rate onto thesecondary accretion shock is so low at late times that it can startto expand after “secondary shock revival”, thus re-establishing anoutflow (bottom right panel of Figure 23).While not indicative of a failure of the explosion, the smallor vanishing outflow surface fraction in the long-time simulationsnonetheless indicates (like the models of Bruenn et al. 2014) thatthe separation of outgoing and infalling mass shells in 2D worksdi ff erently from the usual picture where a high-entropy neutrino-driven wind with an approximately spherical flow geometry even-tually develops. The polar outflows can be viewed as a confinedwind driven jointly by neutrino heating and acoustic waves, butthey never cover the entire sphere, and because of their flow geom-etry and the strong activity of acoustic waves, the outflow dynamicsis completely di ff erent from spherical winds driven purely by neu-trino heating.From the foregoing, it is clear that the inhibition of theneutrino-driven wind in 2D is probably largely artificial; and weonly mention this peculiarity for that reason. As discussed in Sec-tion 5.2.1, the presence of a larger e ff ective eddy viscosity in 3Dcould terminate accretion earlier than in 2D (where supersonicallyinfalling matter is hardly decelerated by lateral momentum trans-fer), or at least decelarate infalling matter su ffi ciently to be sweptalong by an incipient spherical wind after a few seconds. Moreover,our models likely underestimate the di ff usive neutrino luminosity c (cid:13)000 , 000–000 B. M¨uller
20 100 1000r [km]-5-4-3-2-1012 r ∫ δ P δ v r d Ω [ e r g s - ] 〈 ( v r − 〈 v r 〉 ) 〉 [ c m s - ] gain radius Figure 20.
Top: Radial velocity dispersion (cid:104) ( v r − (cid:104) v r (cid:105) ) (cid:105) in 2D and 3D ata post-bounce time of 400 ms. Bottom: Radial profiles of the “acoustic”energy flux r (cid:82) δ P δ v r d Ω in 2D and 3D at a post-bounce time of 400 ms.The curves show temporal averages over several time steps. from the neutron star core and hence the neutrino heating at latetimes because we ignore the e ff ect of nucleon correlations (Bur-rows & Sawyer 1998, 1999; Reddy et al. 1999), which shorten theproto-neutron star cooling time-scale considerably (H¨udepohl et al.2009). It is conceivable that the concomitant increase of the windmass loss rate could still lead to a volume-filling outflow for morerealistic neutrino opacities after a few seconds even in 2D. Based on a successful supernova simulation of a 9 . M (cid:12) star, Mel-son, Janka & Marek (2015) recently suggested that the more ef-ficient braking of the accretion downflows can be responsible forslightly higher explosion energies in 3D because the less violentimpact of the downflows on the neutron star surface lead to reducedcooling. In our comparison of models s11.2 2Da and s11.2 3D weobserve some of the same symptoms noted by these authors, i.e. tur-bulent braking of the downflows and a reduced cooling rate ˙ Q cool atlate times (Figure 13). Obviously, this raises the question whetherthe mechanism proposed by Melson, Janka & Marek (2015) alsooperates in our 3D simulation, and how the physical processes wediscussed so far in Sections 5.2.1 to 5.2.4 are related to it. Unfor-tunately, a comparison with Melson, Janka & Marek (2015) is notstraightforward. While they found a sizable increase of the explo-sion energy of 10% in 3D compared to 2D, their explanation in-volved relatively tiny di ff erences in some quantities (e.g. the gain s u rf ace fr ac ti on o f ou t f l o w s a t r = m s11.2_2Das11.2_3D0 5 10time after bounce [s]00.20.40.60.8 s u rf ace fr ac ti on o f ou t f l o w s a t r = m s11.0_2Ds11.2_2Dbs11.4_2Ds11.6_2D Figure 21.
Top: Surface fraction occupied by the outflows in modelss11.2 2Da (black) and s11.2 3D (red) at a radius of 400 km. The surfacefraction is relatively stable with some fluctuations around 0 . radius and the temperature profiles in 2D and 3D) that cannot beconfidently diagnosed in simulations like ours where the 2D and3D models start to deviate from each other already shortly afterbounce once prompt convection develops (which was not simu-lated by Melson, Janka & Marek 2015). Nonetheless, there is suf-ficient evidence that we observe some rather di ff erent phenomenathan Melson, Janka & Marek (2015).Essentially, the mechanism proposed by Melson, Janka &Marek (2015) involves a recession of the gain radius in 3D com-pared to 2D due to reduced cooling to eject slightly more mate-rial in the explosion. Our simulations agree with Melson, Janka &Marek (2015) in showing a smaller volume-integrated cooling ratein 3D in the long term as accretion slowly subsides (Figure 13).However, we do no find a faster recession of the gain radius in3D during the explosion phase (middle panel of Figure 14), and thesituation is ambiguous for the temperature at the gain radius T gain (bottom panel of Figure 14). In 3D, the temperature T gain stagnatesand falls below the 2D value around 250 ms at a time when theexplosion is already considerably more vigorous in 3D than in 2D.We believe that the stagnation of T gain is more indicative of the slower recession of the gain radius rather than of a higher coolinge ffi ciency: The higher values of the Bernoulli integral and the totalenergy of the downflows at the gain radius in 3D (Figure 17 and18) c (cid:13) , 000–000 ynamics of Neutrino-Driven Supernova Explosions Figure 19.
Left: Entropy s in k b / nucleon in the vicinity of the proto-neutron star in 2D (left half of panel) and 3D (right half of panel) at a post-bounce time of400 ms (identical to Figure 16 except for the zoom level). Right: Heating / cooling rate in MeV / nucleon in 2D (left half of panel) and 3D (right half of panel).Iso-velocity contours for a radial velocity of v r = − cm s − are shown to indicate the location of the accretion downflows. Note that much of the neutrinoheating occurs in the downflows and the confined high-entropy bubble in the equatorial region and hence does not drive an outflow. imply that there is actually more energy per unit mass availablethat can be radiated away in neutrinos as the accreted matter settlesdown in the cooling region.Instead, the faster decline of the accretion rate onto the proto-neutron star in 3D is the dominant factor behind the lower coolingrate, making the lower cooling rate a symptom rather than a cause of the more vigorous explosion. Detailed comparisons would berequired to check whether this true for the 9 . M (cid:12) model of Mel-son, Janka & Marek (2015) as well. Since outflow constriction isunlikely to happen for a model with robust shock expansion, the2D /
3D di ff erences found by Melson, Janka & Marek (2015) as wellas in the parameterize simulations of Handy, Plewa & Odrzywołek(2014) are probably most closely related to the di ff erent outflow ef-ficiency in 2D and 3D, i.e. a more e ffi cient “rerouting” of freshlyaccreted matter into outflows. This tallies with their finding of asmaller surface filling factor of the downflows in 3D, which im-plies that a smaller fraction of the neutrino heating is wasted inregions where it cannot directly power an outflow. It also accountsfor reduced mass accretion into the gain region and hence a reces-sion of the gain radius in mass coordinate. While this mechanismis similar to the one discussed in Section 5.2.3 for our models, thee ff ect is apparently smaller in the simulations of Melson, Janka &Marek (2015) because the accretion subsides fast enough to avoidthe formation of secondary shocks and confined high-entropy bub-bles in 2D, which can reduce the outflow e ffi ciency by a factor of ∼ /
3D di ff erences in the simulations of Mel-son, Janka & Marek (2015). While they take reduced convectiveovershooting in 3D as an indication for reduced wave excitation,the e ff ect probably plays a minor role in their case. The rela-tively small average speeds of the downflows at the gain radius and( ∼ cm s − compared to ∼ cm s − in our model) are boundto make the excitation of g -modes rather ine ffi cient and thus rulethem out as a major energy drain on the gain region in 2D. This is also suggested by the fact they find similar internal energies (andhence binding energies) in the outer regions of the cooling layer in2D and 3D despite the stronger recession of the gain radius in 3D,which is quite di ff erent from what we discussed in Sections 5.2.1and 5.2.2. We have presented a successful 3D GR simulation of the explosionof an 11 . M (cid:12) star using the FMT multi-group transport schemeof M¨uller & Janka (2015). The model has been evolved to almost1 s after bounce, and has reached a diagnostic explosion energy of1 . × erg at that point, which is still growing by the end of thesimulation. The baryonic neutron star mass M by at the end of thesimulations has reached 1 . M (cid:12) and estimates of the final neutronstar mass yield M by ≈ . . . . . M (cid:12) and a gravitational mass notexceeding 1 . M (cid:12) , which is compatible with the measured neutronstar mass distribution (Schwab, Podsiadlowski & Rappaport 2010).The fact that we obtain an explosion for this progenitor with a rela-tively accurate multi-group transport scheme further illustrates thateven the non-exploding state-of-the-art models (Hanke et al. 2013;Tamborra et al. 2014a,b) with the best available neutrino transportand microphysics are apparently very close to shock revival, some-thing which is also suggested by the recent successful 3D explosionmodels of the Garching (Melson, Janka & Marek 2015; Melsonet al. 2015) and Oakridge groups (Lentz et al. 2015).A comparison of the explosion dynamics after shock revivalwith 2D long-time simulations of di ff erent progenitors with ZAMSmasses between 11 . M (cid:12) and 11 . M (cid:12) revealed a faster and morestable growth of the explosion energy in 3D compared to 2D. Be-cause accretion downflows and neutrino-driven outflows coexistover several seconds in 2D, the explosion energy in the 2D mod-els can still reach values of up to 2 × erg, but this comes at theexpense of high neutron star masses ( M by (cid:38) . M (cid:12) ) that are likely c (cid:13)000
3D di ff erences in the simulations of Mel-son, Janka & Marek (2015). While they take reduced convectiveovershooting in 3D as an indication for reduced wave excitation,the e ff ect probably plays a minor role in their case. The rela-tively small average speeds of the downflows at the gain radius and( ∼ cm s − compared to ∼ cm s − in our model) are boundto make the excitation of g -modes rather ine ffi cient and thus rulethem out as a major energy drain on the gain region in 2D. This is also suggested by the fact they find similar internal energies (andhence binding energies) in the outer regions of the cooling layer in2D and 3D despite the stronger recession of the gain radius in 3D,which is quite di ff erent from what we discussed in Sections 5.2.1and 5.2.2. We have presented a successful 3D GR simulation of the explosionof an 11 . M (cid:12) star using the FMT multi-group transport schemeof M¨uller & Janka (2015). The model has been evolved to almost1 s after bounce, and has reached a diagnostic explosion energy of1 . × erg at that point, which is still growing by the end of thesimulation. The baryonic neutron star mass M by at the end of thesimulations has reached 1 . M (cid:12) and estimates of the final neutronstar mass yield M by ≈ . . . . . M (cid:12) and a gravitational mass notexceeding 1 . M (cid:12) , which is compatible with the measured neutronstar mass distribution (Schwab, Podsiadlowski & Rappaport 2010).The fact that we obtain an explosion for this progenitor with a rela-tively accurate multi-group transport scheme further illustrates thateven the non-exploding state-of-the-art models (Hanke et al. 2013;Tamborra et al. 2014a,b) with the best available neutrino transportand microphysics are apparently very close to shock revival, some-thing which is also suggested by the recent successful 3D explosionmodels of the Garching (Melson, Janka & Marek 2015; Melsonet al. 2015) and Oakridge groups (Lentz et al. 2015).A comparison of the explosion dynamics after shock revivalwith 2D long-time simulations of di ff erent progenitors with ZAMSmasses between 11 . M (cid:12) and 11 . M (cid:12) revealed a faster and morestable growth of the explosion energy in 3D compared to 2D. Be-cause accretion downflows and neutrino-driven outflows coexistover several seconds in 2D, the explosion energy in the 2D mod-els can still reach values of up to 2 × erg, but this comes at theexpense of high neutron star masses ( M by (cid:38) . M (cid:12) ) that are likely c (cid:13)000 , 000–000 B. M¨uller
Figure 22.
Constriction and partial shredding of an outflow in model s11.2 2Da, shown by snapshots of the entropy at post-bounce times of 592 ms, 611 ms,628 ms, and 655 ms. A downflow (cyan) originating from a Rayleigh-Taylor plume of cold matter penetrates the hot-neutrino heated bubble in the northernhemispheres (top left), constricts the neutrino-heated bubble to a tenuous outflow as it approaches the axis (top right), and eventually a considerable amount ofcold material is mixed into the outflow (bottom left). While the ejection of matter continues (bottom right), the mixing event lowers the average total energyper unit mass in the ejecta. incompatible with the observed neutron star mass distribution. Adetailed comparison of the 2D and 3D models unearthed severalphysical mechanisms responsible for the more robust rise of theexplosion energy in 3D and the slower growth of the proto-neutronstar mass, which is a symptom of the faster subsidence of accre-tion. The specific e ff ects that we find are summarized below, andwe also provide a schematic visualization of the di ff erent flow ge-ometry and the energy budget between the downflows, the outflowsand the gain region in Figure 24 to aid the reader’s understanding:1. In 2D, the interfaces between the accretion funnels and theneutrino-heated bubbles tend to become laminar after shock re-vival, while they are corrugated by the Kelvin-Helmholtz instabilityin 3D. We ascribe the di ff erent behavior to the suppression of thepurely two-dimensional modes of the Kelvin-Helmholtz instabil-ity in the supersonic regime (Gerwin 1968). The e ff ect is thus dis-tinct from the inverse turbulent energy cascade in 2D (Kraichnan1967), which has been invoked as an explanation for the di ff erentbehaviour of 2D and 3D models prior to shock revival, since the di ff erent energy cascade in 2D and 3D is not related to the Machnumber of the flow.2. As a consequence, the e ff ective eddy viscosity and di ff usiv-ity between the downflows and outflows is larger in 3D than in 2Dduring most phases, i.e. there is more exchange of energy and mo-mentum between the outflows and downflows. On the one hand,this implies that the outflows contribute only ∼ / nucleonto the explosion energy in 3D, as some of the net total (i.e ther-mal + kinetic + potential) energy gained from nucleon recombinationof ∼ . / nucleon is lost to the downflows by turbulent dif-fusion. On the other hand, the turbulence e ff ectively “brakes” thedownflows, and they arrive at the gain radius with smaller velocitiesbut higher total energy per unit mass than in 2D.3. The higher impact velocities of the downflows and the for-mation of secondary accretion shocks at small radii in 2D lead to amore e ffi cient excitation of g -modes and acoustic waves at the gainradius that transport energy into deeper regions of the cooling layerand into the neutrino-heated ejecta respectively. Our analysis sug- c (cid:13) , 000–000 ynamics of Neutrino-Driven Supernova Explosions Figure 23.
Snapshots of the radial velocity in units of 10 cm (left half of panels, lower color bar) and the specific entropy s in k b / nucleon (right half ofpanels, upper color bar) depicting the separation of the high-entropy outflow from the gain region in model s11.6 2D and the re-establishment of an outflowafter “secondary shock revival” in model s11.0 2D. Red isovelocity contours are used to separate outward-moving matter with radial velocities larger than10 cm s − from infalling matter. At 4 . . . ff earlier. By the end ofthe simulation at 8 . . M (cid:12) in the downflows are still falling toward the proto-neutron star. gests that the energy loss from the gain region by wave excitationbecomes comparable to the volume-integrated neutrino heating rateat late times, and by increasing the absolute value of the bindingenergy | e tot | at the gain radius significantly reduce the mass outflowrate that can be sustained by neutrino heating. In 3D, the turbulentenergy flux into the gain region is small, and the binding energyat the gain radius is smaller by factor of (cid:38) . / nucleon.4. In addition, the outflow e ffi ciency η out = ˙ M out / ( ˙ Q heat / | e gain | ), is also higher in 3D ( η out ∼ η out (cid:46) . ffi ciency in 2D stems from the large sur-face fraction occupied by fast downflows and “confined bubbles”between downflows whose expansion is inhibited by the formationof secondary accretion shocks.5. Episodic mixing between the outflows and downflows stilloccurs in 2D, e.g. by the formation of new downflows as a resultof the Rayleigh-Taylor between the cold shocked material and theneutrino-heated high entropy bubbles. While mixing only occurssporadically in 2D, the consequences of these mixing events are c (cid:13)000
Snapshots of the radial velocity in units of 10 cm (left half of panels, lower color bar) and the specific entropy s in k b / nucleon (right half ofpanels, upper color bar) depicting the separation of the high-entropy outflow from the gain region in model s11.6 2D and the re-establishment of an outflowafter “secondary shock revival” in model s11.0 2D. Red isovelocity contours are used to separate outward-moving matter with radial velocities larger than10 cm s − from infalling matter. At 4 . . . ff earlier. By the end ofthe simulation at 8 . . M (cid:12) in the downflows are still falling toward the proto-neutron star. gests that the energy loss from the gain region by wave excitationbecomes comparable to the volume-integrated neutrino heating rateat late times, and by increasing the absolute value of the bindingenergy | e tot | at the gain radius significantly reduce the mass outflowrate that can be sustained by neutrino heating. In 3D, the turbulentenergy flux into the gain region is small, and the binding energyat the gain radius is smaller by factor of (cid:38) . / nucleon.4. In addition, the outflow e ffi ciency η out = ˙ M out / ( ˙ Q heat / | e gain | ), is also higher in 3D ( η out ∼ η out (cid:46) . ffi ciency in 2D stems from the large sur-face fraction occupied by fast downflows and “confined bubbles”between downflows whose expansion is inhibited by the formationof secondary accretion shocks.5. Episodic mixing between the outflows and downflows stilloccurs in 2D, e.g. by the formation of new downflows as a resultof the Rayleigh-Taylor between the cold shocked material and theneutrino-heated high entropy bubbles. While mixing only occurssporadically in 2D, the consequences of these mixing events are c (cid:13)000 , 000–000 B. M¨uller
Figure 24.
Sketch of the di ff erent energy budget between the outflows (yellow), the downflows (red) and the cooling region (blue) in 2D (left) and 3D (right).In 3D, turbulent eddy di ff usivity leads to a persistent, small energy flux (short solid arrow) from the neutrino-heated outflows to the downflows, whereasmixing between the outflows and downflows only occurs episodically in 2D, but has more dramatic consequences in this case because it leads to large-scalemixing of cold matter into the outflows (long dotted arrow). There is a considerable energy transfer from the gain region to the cooling region due to waveexcitation in 2D and an indirect transfer of energy from the downflows to the gain region by the excitation of acoustic waves in 2D (which can lead to largerasymptotic energies per unit mass in the ejecta); this is absent in 3D. Moreover, a considerable amount of neutrino heating (red arrows) is wasted in 2D becauseit is deposited in confined bubbles, whereas almost the entire neutrino heating in 3D is used to lift matter out of the gravitational well. more catastrophic than in 3D. Not only do they slow down the riseof the explosion energy by mixing cold material into the outflows;the penetration of accretion funnels into the high-entropy bubblescan also lead to the constriction of outflows, sometimes shuttingthem o ff completely and permanently decreasing the outflow sur-face fraction to values of (cid:46) . ff ectscan play a beneficial role in core-collapse supernova explosions af-ter shock revival. However, since our current 3D model has onlybeen evolved to ∼ ff ects duringthe explosion phase can be summarized as follows:1. Longer 3D simulations with higher resolution are necessaryto determine final explosion energies, Nickel masses and neutronstar masses precisely for comparison with observations without re-course to indirect methods. Our estimate for the final baryonic neu-tron star mass of 1 . M (cid:12) for the 3D model of the 11 . M (cid:12) pro-genitor still assumes the accretion of an additional 0 . M (cid:12) solarmasses, which is much more than a visual inspection of Figure 10suggests given the very slow rise of the neutron star mass in 3D. Ifthe turbulent braking of the downflows terminates accretion beforethe post-shock velocity equals the escape velocity as speculatedin Section 5.2.1, the final baryonic and gravitational neutron starmass might be as low as ∼ . M (cid:12) and ∼ . M (cid:12) , respectively. Thiswould indicate that a plausible distribution of neutron star massesspanning the entire range of observed values down to the lower endis within reach of modern multi-D simulations of neutrino-drivensupernovae.2. A more rigorous analysis of the turbulent multi-dimensionalflow in the spirit of Reynolds decomposition would be highly de-sirable in order to further bolster our qualitative interpretation of3D e ff ects in the post-explosion phase. Such quantitative analysismethods have considerably advanced our understanding of the tur- bulent flow during the accretion phase (Murphy & Meakin 2011).After shock revival, the nonstationarity of the flow presents a chal-lenge for such methods, however. Our relatively crude two-streamanalysis based on a separation of the outflows and downflows couldalso be improved in order to account more directly and rigorouslyfor the turbulent exchange of mass, momentum and energy betweenthe two streams, but such an analysis faces a major challenge in theform of the complicated flow geometry.3. Whether and to what extent the positive 3D e ff ects describedin this paper come into play obviously depends on whether shockrevival can be accomplished in 3D in the first place and on the de-lay compared to the 2D case. If there is a significant delay in shockrevival, 3D models may not be able to equalize the “head start” ofthe 2D models at least of relatively powerful explosions where thediagnostic energy shows first signs of levelling o ff after ∼
300 msor less (Bruenn et al. 2014; Pan et al. 2015). Even in this case, themechanism discussed in this paper could nonetheless help to mit-igate the “penalty” incurred by the delay of the explosion in 3Dand allow the models to remain compatible with observational con-straints. Moreover, if accretion lasts longer — as in the 2D simu-lations of M¨uller, Janka & Marek (2012); M¨uller, Janka & Heger(2012); Janka et al. (2012) — the beneficial 3D e ff ects in the phaseafter shock revival may outweigh the “penalty” of delayed shockrevival. Furthermore, it is conceivable that the problem of missingor delayed explosion in 3D may yet be resolved by the inclusion ofbetter, multi-dimensional progenitor models with large-scale initialperturbations that aid shock revival (Couch & Ott 2013; M¨uller &Janka 2015; Couch & Ott 2015), unknown microphysics (Melsonet al. 2015); and strongly SASI-dominated models may even ex-plode easier in 3D (Fern´andez 2015).4. The robustness of the mechanisms described in this paperneeds to be studied further for a wider range of progenitors. The2D simulations (Buras et al. 2006a; M¨uller, Janka & Marek 2012)of the 11 . M (cid:12) progenitor considered here have been particularlynoteworthy examples for suspiciously low explosion energies andlong-lasting accretion. This behavior is due to the specific charac-teristics of progenitors around 11 M (cid:12) , including a relatively smallsilicon core and a very pronounced density jump at the Si / SiO in-terface. These properties result in a small proto-neutron star mass M immediately after shock revival and hence low neutrino energies c (cid:13) , 000–000 ynamics of Neutrino-Driven Supernova Explosions (cp. the scaling of the electron antineutrino energy with M foundby M¨uller & Janka 2014) as well as a small accretion luminosity.Both of these factors contribute to relatively weak neutrino heatingafter shock revival and a small mass outflow rate. The tepid natureof our 2D explosions may thus hinge very much on the peculiarstructure of low-mass supernova progenitors.It therefore remains to be seen whether 3D e ff ects provide a sim-ilarly strong boost for the growth of the explosion energy in otherprogenitors. Whenever 2D models develop the characteristic broaddownflows and secondary shocks indicative of long-lasting accre-tion during a relatively weak explosions, such as the 15 M (cid:12) and27 M (cid:12) models of M¨uller, Janka & Marek (2012) and Janka et al.(2012), the physical mechanisms identified here likely come intoplay eventually. On the other hand, they may play a negligible roleif the volume fraction of the downflows drops very quickly as in themodels of Bruenn et al. (2014) and Lentz et al. (2015) or the pa-rameterized models of Handy, Plewa & Odrzywołek (2014). Since2D supernova simulations of di ff erent groups have not yet con-verged su ffi ciently to decide whether there is a generic problemof “weak explosions” in 2D, it is impossible to judge the genericcharacter of our findings. By the same token, however, it cannotbe ruled out that 2D models should be generically underenergeticand overestimate the amount of accretion after shock revival due tothe mechanisms we identified, and that realistic explosion enegiesand remnant masses will only be obtained in 3D. If so, prematurelyconfronting the 2D models with the observational constraints couldlead to wrong conclusions.Thus, more work is necessary to substantiate the intriguing per-spective that 3D e ff ects could help to achieve agreement betweensimulations of neutrino-driven supernovae and observational con-straints such as explosion energies and neutron star masses. Farfrom o ff ering a complete solution due to the limitations of compu-tational resources that have always plagued supernova theory, ourpresent study can only take a first step in this direction and adum-brate some of the physical mechanisms that could help to boost3D explosions after shock revival. Nonetheless, even our currentresults already serve an antidote against undue pessimism after ini-tial setbacks in 3D multi-group neutrino hydrodynamics simula-tions. Along with the recent successful explosions in first-principlemodels, the identification of other beneficial e ff ects of the third di-mension on the explosion threshold and energetics, and plausibleideas for solving the problem of missing explosions with the helpof multi-D progenitor models and / or non-standard microphysics,they are another piece that fits well into the overall puzzle, sug-gesting that a solution for the supernova problem is slowly takingshape. ACKNOWLEDGEMENTS
We acknowledge fruitful exchange with H. Andressen, A. Burrows,Th. Foglizzo, A. Heger, W.R. Hix, H.-Th. Janka, Y. Levin, T. Plewa,and T. Waters. The author has been supported by the AustralianResearch Council through a Discovery Early Career ResearcherAward (grant DE150101145) and by the Alexander von HumboldtFoundation through a Feodor Lynen fellowship. The computationswere performed on
Raijin at the NCI National Facility (project fh6)using computer time contingents obtained through NCMAS, AS-TAC, and a Monash LIEF top-up grant, on the Monash eGrid Clus-ter, and on the IBM iDataPlex system hydra at the Rechenzentrumof the Max-Planck Society (RZG).
APPENDIX A: DEFINITION OF ENERGIES ANDENERGY FLUXES IN GENERAL RELATIVITY
While the computation of mass fluxes and spherically-averagedprofiles of hydrodynamic quantities can be readily generalized tothe relativistic case just by including the correct three-volume ele-ment, the definition of energies and energy fluxes in the relativisticcase is less straightforward, and is therefore briefly expounded inthis Appendix. We use geometric units ( G = c =
1) throughout thissection.After adopting a 3 + τ as ∂ √ γτ∂ t + ∂ √− g (cid:16) τ ˆ v i + Pv i (cid:17) ∂ x i = α √− g (cid:32) T µ ∂ ln α∂ x µ − T µν Γ µν (cid:33) . (A1)Here, τ is defined in terms of the baryonic mass density ρ in thefluid frame, the Lorentz factor W , the internal energy density (cid:15) (in-cluding all rest-mass contributions), the pressure P , and the rela-tivistic enthalpy h = + (cid:15) + P /ρ as τ = ρ hW − P − ρ W = ρ (1 + (cid:15) + P /ρ ) W − P − ρ W . (A2)Furthermore, γ and g are the determinants of the three- and four-metric, respectively, and the advection term contains ˆ v i = v i − β i /α instead of the Eulerian three-velocity v i . By pushing the lapse func-tion α into the temporal and spatial derivatives, it is possible to for-mulate a strict conservation law (without a source term) analogousto in the limit of a stationary space-time with a zero shift vector(cp. Equation A35 M¨uller, Janka & Dimmelmeier 2010, where theright-hand side reduces to zero in this limit): ∂∂ t (cid:2) √ γα ( τ + D ) − √ γ D (cid:3) + ∂∂ x i (cid:104) √− g (cid:16) ατ ˆ v i + α D ˆ v i − D ˆ v i + α Pv i (cid:17)(cid:105) = . (A3)Here, we have introduced the baryonic mass density in the Eulerianframe D = ρ W to simplify the equation.This suggests that in the limit of a vanishing shift vector and astationary metric, the Newtonian expression e tot = (cid:15) + v / + Φ forthe total energy per unit mass (including rest-mass contributions)can be generalized to e tot , rm = ατ D + ( α − , (A4)and the role of the total Newtonian enthalpy in the flux is taken by h tot , rm = ατ/ D + ( α − + α P / D . (A5)It is noteworthy that the internal energy and rest-mass contributions(which enter the equations through τ ) always appears in conjunc-tion with factors W and α . Strictly speaking, it is therefore no longerpossible to formulate the energy equation in general relativity with-out including rest-mass contributions in the conserved quantitiesand the fluxes by pushing them into a nuclear source term instead(as least not in a simple form). Computing fluxes and total ener-gies excluding rest-mass contributions is therefore somewhat lessmeaningful in general relativity. However, since we have α ≈ v (cid:46) . c in the gain region, the higher-order relativistic correctionsare small enough to be neglected, it is still reasonable to computetotal energies and enthalpies without rest-mass contributions by us-ing just the thermal energy contribution (cid:15) therm to the internal energydensity (cid:15) instead of (cid:15) = (cid:15) therm + (cid:15) rm .For the other quantities, considered in this paper, the gener-alization is trivial. Mass fluxes through the surface of a sphere or c (cid:13)000
1) throughout thissection.After adopting a 3 + τ as ∂ √ γτ∂ t + ∂ √− g (cid:16) τ ˆ v i + Pv i (cid:17) ∂ x i = α √− g (cid:32) T µ ∂ ln α∂ x µ − T µν Γ µν (cid:33) . (A1)Here, τ is defined in terms of the baryonic mass density ρ in thefluid frame, the Lorentz factor W , the internal energy density (cid:15) (in-cluding all rest-mass contributions), the pressure P , and the rela-tivistic enthalpy h = + (cid:15) + P /ρ as τ = ρ hW − P − ρ W = ρ (1 + (cid:15) + P /ρ ) W − P − ρ W . (A2)Furthermore, γ and g are the determinants of the three- and four-metric, respectively, and the advection term contains ˆ v i = v i − β i /α instead of the Eulerian three-velocity v i . By pushing the lapse func-tion α into the temporal and spatial derivatives, it is possible to for-mulate a strict conservation law (without a source term) analogousto in the limit of a stationary space-time with a zero shift vector(cp. Equation A35 M¨uller, Janka & Dimmelmeier 2010, where theright-hand side reduces to zero in this limit): ∂∂ t (cid:2) √ γα ( τ + D ) − √ γ D (cid:3) + ∂∂ x i (cid:104) √− g (cid:16) ατ ˆ v i + α D ˆ v i − D ˆ v i + α Pv i (cid:17)(cid:105) = . (A3)Here, we have introduced the baryonic mass density in the Eulerianframe D = ρ W to simplify the equation.This suggests that in the limit of a vanishing shift vector and astationary metric, the Newtonian expression e tot = (cid:15) + v / + Φ forthe total energy per unit mass (including rest-mass contributions)can be generalized to e tot , rm = ατ D + ( α − , (A4)and the role of the total Newtonian enthalpy in the flux is taken by h tot , rm = ατ/ D + ( α − + α P / D . (A5)It is noteworthy that the internal energy and rest-mass contributions(which enter the equations through τ ) always appears in conjunc-tion with factors W and α . Strictly speaking, it is therefore no longerpossible to formulate the energy equation in general relativity with-out including rest-mass contributions in the conserved quantitiesand the fluxes by pushing them into a nuclear source term instead(as least not in a simple form). Computing fluxes and total ener-gies excluding rest-mass contributions is therefore somewhat lessmeaningful in general relativity. However, since we have α ≈ v (cid:46) . c in the gain region, the higher-order relativistic correctionsare small enough to be neglected, it is still reasonable to computetotal energies and enthalpies without rest-mass contributions by us-ing just the thermal energy contribution (cid:15) therm to the internal energydensity (cid:15) instead of (cid:15) = (cid:15) therm + (cid:15) rm .For the other quantities, considered in this paper, the gener-alization is trivial. Mass fluxes through the surface of a sphere or c (cid:13)000 , 000–000 B. M¨uller parts of it are computed as˙ M = (cid:90) α Dv r r φ d Ω , (A6)where v r is the radial velocity component measured in the orthonor-malised Eulerian frame and φ is the conformal factor in the xCFCmetric; and the computation of energy / enthalpy fluxes works anal-ogously. The density-weighted spherical average ¯ X of a quantity X is computed as ¯ X = (cid:82) DX φ d Ω (cid:82) D φ d Ω . (A7) References
Abdikamalov E. et al., 2014, ArXiv e-prints, 1409.7078Arnett D., 1994, ApJ, 427, 932Arnett W. D., Meakin C., 2011, ApJ, 733, 78Asida S. M., Arnett D., 2000, ApJ, 545, 435Balsa T. F., Goldstein M. E., 1990, Journal of Fluid Mechanics,216, 585Banyuls F., Font J. A., Ibanez J. M. A., Mart´ı J. M. A., MirallesJ. A., 1997, ApJ, 476, 221Bazan G., Arnett D., 1994, ApJL, 433, L41Bazan G., Arnett D., 1998, ApJ, 496, 316Blondin J. M., Mezzacappa A., DeMarino C., 2003, ApJ, 584, 971Blumen W., 1970, Journal of Fluid Mechanics, 40, 769Blumen W., Drazin P. G., Billings D. F., 1975, Journal of FluidMechanics, 71, 305Boyd J. P., 2001, Chebyshev and Fourier Spectral Methods, 2nded. Dover Publications, MineolaBruenn S. W., 1985, ApJS, 58, 771Bruenn S. W. et al., 2014, ArXiv e-prints, 1409.5779Bruenn S. W. et al., 2013, ApJL, 767, L6Buras R., Janka H.-T., Rampp M., Kifonidis K., 2006a, A&A,457, 281Buras R., Rampp M., Janka H.-T., Kifonidis K., 2006b, A&A,447, 1049Burrows A., 2013, Reviews of Modern Physics, 85, 245Burrows A., Fryxell B. A., 1992, Science, 258, 430Burrows A., Goshy J., 1993, ApJL, 416, L75 + Burrows A., Hayes J., Fryxell B. A., 1995, ApJ, 450, 830Burrows A., Livne E., Dessart L., Ott C. D., Murphy J., 2006,ApJ, 640, 878Burrows A., Livne E., Dessart L., Ott C. D., Murphy J., 2007,ApJ, 655, 416Burrows A., Sawyer R. F., 1998, Phys. Rev. C, 58, 554Burrows A., Sawyer R. F., 1999, Phys. Rev. C, 59, 510Cerd´a-Dur´an, 2009, Mesh coarsening in relativistic hydro-dynamics,
Choudhury S. R., Lovelace R. V. E., 1984, ApJ, 283, 331Chugai N. N., Utrobin V. P., 2014, Astronomy Letters, 40, 291Colella P., Woodward P. R., 1984,
J. Comp. Phys. , 54, 174Cordero-Carri´on I., Cerd´a-Dur´an P., Dimmelmeier H., JaramilloJ. L., Novak J., Gourgoulhon E., 2009, Phys. Rev. D, 79, 024017Couch S. M., 2013a, ApJ, 775, 35Couch S. M., 2013b, ApJ, 765, 29Couch S. M., Chatzopoulos E., Arnett W. D., Timmes F. X., 2015,ArXiv e-prints, 1503.02199Couch S. M., O’Connor E. P., 2014, ApJ, 785, 123 Couch S. M., Ott C. D., 2013, ApJL, 778, L7Couch S. M., Ott C. D., 2015, ApJ, 799, 5Dimmelmeier H., Font J. A., M¨uller E., 2002, A&A, 388, 917Dolence J. C., Burrows A., Murphy J. W., Nordhaus J., 2013, ApJ,765, 110Drazin P. G., Davey A., 1977, Journal of Fluid Mechanics, 82, 255Fern´andez R., 2015, ArXiv e-prints, 1504.07996Foglizzo T., Galletti P., Scheck L., Janka H.-T., 2007, ApJ, 654,1006Fragile P. C., Lindner C. C., Anninos P., Salmonson J. D., 2009,ApJ, 691, 482Gerwin R. A., 1968, Reviews of Modern Physics, 40, 652Goldreich P., Kumar P., 1990, ApJ, 363, 694Guilet J., Foglizzo T., 2012, MNRAS, 421, 546Handy T., Plewa T., Odrzywołek A., 2014, ApJ, 783, 125Hanke F., 2014, PhD thesis, Technische Universti¨at M¨unchenHanke F., Marek A., M¨uller B., Janka H.-T., 2012, ApJ, 755, 138Hanke F., M¨uller B., Wongwathanarat A., Marek A., Janka H.-T.,2013, ApJ, 770, 66Herant M., Benz W., Colgate S., 1992, ApJ, 395, 642Herant M., Benz W., Hix W. R., Fryer C. L., Colgate S. A., 1994,ApJ, 435, 339H¨udepohl L., 2014, PhD thesis, Technische Universti¨at M¨unchenH¨udepohl L., M¨uller B., Janka H., Marek A., Ra ff elt G. G., 2009,Phys. Rev. Lett.Janka H.-T., 2012, Annual Review of Nuclear and Particle Sci-ence, 62, 407Janka H.-T., Hanke F., H¨udepohl L., Marek A., M¨uller B., Ober-gaulinger M., 2012, Progress of Theoretical and ExperimentalPhysics, 2012, 010000Janka H.-T., M¨uller E., 1996, A&A, 306, 167Kageyama A., Sato T., 2004, Geochemistry, Geophysics, Geosys-tems, 5, n / a, q09005Kifonidis K., Plewa T., Janka H.-T., M¨uller E., 2003, A&A, 408,621Klimishin I. A., Gnatyk B. I., 1981, Astrophysics, 17, 306Koldoba A. V., Romanova M. M., Ustyugova G. V., LovelaceR. V. E., 2002, ApJL, 576, L53Kompaneets A. S., 1960, Soviet Physics Doklady, 5, 46Koo B.-C., McKee C. F., 1990, ApJ, 354, 513Kraichnan R. H., 1967, Physics of Fluids, 10, 1417Kuhlen M., Woosley W. E., Glatzmaier G. A., 2003, in Astro-nomical Society of the Pacific Conference Series, Vol. 293, 3DStellar Evolution, Turcotte S., Keller S. C., Cavallo R. M., eds.,p. 147Laming J. M., 2007, ApJ, 659, 1449Lattimer J. M., Swesty F. D., 1991, Nucl. Phys. A , 535, 331Laumbach D. D., Probstein R. F., 1969, Journal of Fluid Mechan-ics, 35, 53Lecoanet D., Quataert E., 2013, MNRAS, 430, 2363Lentz E. J. et al., 2015, ArXiv e-prints, 1505.05110Lentz E. J., Mezzacappa A., Bronson Messer O. E., Hix W. R.,Bruenn S. W., 2012a, ApJ, 760, 94Lentz E. J., Mezzacappa A., Bronson Messer O. E., Liebend¨orferM., Hix W. R., Bruenn S. W., 2012b, ApJ, 747, 73Malamud G. et al., 2013, High Energy Density Physics, 9, 672Marek A., Janka H., 2009, ApJ, 694, 664Marek A., Janka H., M¨uller E., 2009, A&A, 496, 475Matzner C. D., McKee C. F., 1999, ApJ, 510, 379Meakin C. A., Arnett D., 2006, ApJL, 637, L53Meakin C. A., Arnett D., 2007a, ApJ, 665, 690Meakin C. A., Arnett D., 2007b, ApJ, 667, 448 c (cid:13) , 000–000 ynamics of Neutrino-Driven Supernova Explosions Melson T., Janka H.-T., Bollig R., Hanke F., Marek A., MuellerB., 2015, ArXiv e-prints, 1504.07631Melson T., Janka H.-T., Marek A., 2015, ApJL, 801, L24Mignone A., Bodo G., 2005, MNRAS, 364, 126M¨uller B., Janka H., Dimmelmeier H., 2010, ApJS, 189, 104M¨uller B., Janka H.-T., 2014, ApJ, 788, 82M¨uller B., Janka H.-T., 2015, MNRAS, 448, 2141M¨uller B., Janka H.-T., Heger A., 2012, ApJ, 761, 72M¨uller B., Janka H.-T., Marek A., 2012, ApJ, 756, 84M¨uller B., Janka H.-T., Marek A., 2013, ApJ, 766, 43M¨uller E., Janka H.-T., 1997, A&A, 317, 140Murphy J. W., Burrows A., 2008, ApJ, 688, 1159Murphy J. W., Dolence J. C., Burrows A., 2013, ApJ, 771, 52Murphy J. W., Meakin C., 2011, ApJ, 742, 74Murphy J. W., Ott C. D., Burrows A., 2009, ApJ, 707, 1173Nakamura K., Kuroda T., Takiwaki T., Kotake K., 2014, ApJ, 793,45Obergaulinger M., Janka H.-T., Aloy M. A., 2014, MNRAS, 445,3169Pan K.-C., Liebend¨orfer M., Hempel M., Thielemann F.-K., 2015,ArXiv e-prints, 1505.02513Pejcha O., Prieto J. L., 2015, ApJ, 799, 215Perego A., Hempel M., Fr¨ohlich C., Ebinger K., Eichler M.,Casanova J., Liebend¨orfer M., Thielemann F.-K., 2015, ApJ,806, 275Poznanski D., 2013, MNRAS, 436, 3224Rampp M., Janka H.-T., 2002, A&A, 396, 361Reddy S., Prakash M., Lattimer J. M., Pons J. A., 1999,Phys. Rev. C, 59, 2888Ronchi C., Iacono R., Paolucci P. S., 1996, Journal of Computa-tional Physics, 124, 93Scheck L., Kifonidis K., Janka H.-T., M¨uller E., 2006, A&A, 457,963Schwab J., Podsiadlowski P., Rappaport S., 2010, ApJ, 719, 722Sedov L. I., 1959, Similarity and Dimensional Methods in Me-chanics. Academic Press, New YorkSuwa Y., Kotake K., Takiwaki T., Whitehouse S. C., Liebend¨orferM., Sato K., 2010, PASJ, 62, L49 + Suwa Y., Takiwaki T., Kotake K., Fischer T., Liebend¨orfer M.,Sato K., 2013, ApJ, 764, 99Takiwaki T., Kotake K., Suwa Y., 2014, ApJ, 786, 83Tamborra I., Hanke F., Janka H.-T., M¨uller B., Ra ff elt G. G.,Marek A., 2014a, ApJ, 792, 96Tamborra I., Ra ff elt G., Hanke F., Janka H.-T., M¨uller B., 2014b,Phys. Rev. D, 90, 045032Tanaka M. et al., 2009, ApJ, 692, 1131Thompson C., 2000, ApJ, 534, 915Timmes F. X., Woosley S. E., Weaver T. A., 1996, ApJ, 457, 834Utrobin V. P., Chugai N. N., 2011, A&A, 532, A100van Leer B., 1977, Journal of Computational Physics, 23, 263Weinberg N. N., Quataert E., 2008, MNRAS, 387, L64Wilson J. R., Mayle R. W., 1988, Phys. Rep., 163, 63Wongwathanarat A., Hammer N. J., M¨uller E., 2010, A&A, 514,A48Wongwathanarat A., M¨uller E., Janka H.-T., 2015, A&A, 577,A48Woosley S. E., Heger A., Weaver T. A., 2002, Rev. Mod. Phys. ,74, 1015Woosley S. E., Weaver T. A., 1995, ApJS, 101, 181Yakunin K. N. et al., 2010, Classical and Quantum Gravity, 27,194005Zink B., Schnetter E., Tiglio M., 2008, Phys. Rev. D, 77, 103015 c (cid:13)000