The Status of Multi-Dimensional Core-Collapse Supernova Models

Abstract

Models of neutrino-driven core-collapse supernova explosions have matured considerably in recent years. Explosions of low-mass progenitors can routinely be simulated in 1D, 2D, and 3D. Nucleosynthesis calculations indicate that these supernovae could be contributors of some lighter neutron-rich elements beyond iron. The explosion mechanism of more massive stars remains under investigation, although first 3D models of neutrino-driven explosions employing multi-group neutrino transport have become available. Together with earlier 2D models and more simplified 3D simulations, these have elucidated the interplay between neutrino heating and hydrodynamic instabilities in the post-shock region that is essential for shock revival. However, some physical ingredients may still need to be added/improved before simulations can robustly explain supernova explosions over a wide range of progenitors. Solutions recently suggested in the literature include uncertainties in the neutrino rates, rotation, and seed perturbations from convective shell burning. We review the implications of 3D simulations of shell burning in supernova progenitors for the ‘perturbations-aided neutrino-driven mechanism,’ whose efficacy is illustrated by the first successful multi-group neutrino hydrodynamics simulation of an 18 solar mass progenitor with 3D initial conditions. We conclude with speculations about the impact of 3D effects on the structure of massive stars through convective boundary mixing.

1 INTRODUCTION

The explosions of massive stars as core-collapse supernovae (CCSNe) constitute one of the most outstanding problems in modern astrophysics. This is in no small measure due to the critical role of supernova explosions in the history of the Universe. CCSNe figure prominently in the chemical evolution of galaxies as the dominant producers, e.g., of elements between oxygen and the iron group (Arnett 1996; Woosley, Heger, & Weaver 2002), and supernova feedback is a key ingredient in the modern theory of star formation (Krumholz 2014). The properties of neutron stars and stellar-mass black holes (masses, spins, kicks; Özel et al. 2010, 2012; Kiziltan et al. 2013; Antoniadis et al. 2016; Arzoumanian, Chernoff, & Cordes 2002; Hobbs et al. 2005) cannot be understood without addressing the origin of these compact objects in stellar explosions.

Why (some) massive stars explode is, however, a daunting problem in its own right regardless of the wider implications of supernova explosions: The connection of supernovae of massive stars with the gravitational collapse to a neutron star has been postulated more than 80 yr ago (Baade & Zwicky 1934), and the best-explored mechanism for powering the explosion, the neutrino-driven mechanism, has gone through several stages of ‘moulting’ in the 50 yr after its conception by Colgate & White (1966). Yet, the problem of the supernova explosion mechanism still awaits a definitive solution. The rugged path towards an understanding of the explosion mechanism merely reflects that CCSNe are the epitome of a ‘multi-physics’ problem that combines aspects of stellar structure and evolution, nuclear and neutrino physics, fluid dynamics, kinetic theory, and general relativity. We cannot recapitulate the history of the field here and instead refer the reader to the classical and modern reviews of Bethe (1990), Arnett (1996), Mezzacappa (2005), Kotake, Sato, & Takahashi (2006), Janka et al. (2007), Burrows et al. (2007a), Janka (2012), and Burrows (2013) as starting points.

The longevity of the supernova problem should not be misinterpreted: Despite the occasional detour, supernova theory has made steady progress, particularly so during the last few years, which have seen the emergence of mature—and increasingly successful—multi-dimensional first-principle simulations of the collapse and explosion of massive stars as well as conceptual advances in our understanding of the neutrino-driven explosion mechanism and its interplay with multi-dimensional hydrodynamic instabilities.

1.1 The neutrino-driven explosion mechanism in its modern flavour

Before we review these recent advances, it is apposite to briefly recapitulate the basic idea of the neutrino-driven supernova mechanism in its modern guise. Stars with zero-age main sequence (ZAMS) masses above ${\gtrsim }8\,\text{M}_\odot$ and with a helium core mass ${\lesssim }65\,\text{M}_\odot$ (the lower limit for non-pulsational pair-instability supernovae; Heger & Woosley 2002; Heger et al. 2003) develop iron cores that eventually become subject to gravitational instability and undergo collapse on a free-fall timescale. For low-mass supernova progenitors with highly degenerate iron cores, collapse is triggered by the reduction of the electron degeneracy pressure due to electron captures; for more massive stars with higher core entropy and a strong contribution of radiation pressure, photo-disintegration of heavy nuclei also contributes to gravitational instability. Aside from these ‘iron-core supernovae’, there may also be a route towards core collapse from super-AGB stars with O–Ne–Mg cores (Nomoto 1984, 1987; Poelarends et al. 2008; Jones et al. 2013; Jones, Hirschi, & Nomoto 2014; Doherty et al. 2015), where rapid core contraction is triggered by electron captures on ²⁰Ne and ²⁴Mg; ¹ hence this sub-class is designated as ‘electron-capture supernovae’ (ECSNe).

According to modern shell-model calculations (Langanke & Martínez-Pinedo 2000; Langanke et al. 2003), the electron capture rate on heavy nuclei remains high even during the advanced stages of collapse (Langanke et al. 2003) when the composition of the core is dominated by increasingly neutron-rich and massive nuclei. Further deleptonisation during collapse thus reduces the lepton fraction Y _lep to about 0.3 according to modern simulations (Marek et al. 2005; Sullivan et al. 2016) until neutrino trapping occurs at a density of $\mathord {\sim }10^{12} \, {\rm g} \, {\rm cm}^{-3}$ . As a result, the homologously collapsing inner core shrinks (Yahil 1983), and the shock forms at a small enclosed mass of $\mathord {\sim } 0.5\,\text{M}_\odot$ (Langanke et al. 2003; Hix et al. 2003; Marek et al. 2005) after the core reaches supranuclear densities and rebounds (bounces). Due to photodisintegration of heavy nuclei in the infalling shells into free nucleons as well as rapid deleptonisation in the post-shock region once the shock breaks out of the neutrinosphere, the shock stalls a few milliseconds after bounce, i.e., it turns into an accretion shock with negative radial velocity downstream of the shock. Aided by a continuous reduction of the mass accretion rate onto the young proto-neutron star, the stalled accretion shock still propagates outward for $\mathord {\sim } 70 \, {\rm ms}$ , however, and reaches a typical peak radius of $\mathord {\sim } 150 \, {\rm km}$ before it starts to recede again.

The point of maximum shock expansion is roughly coincident with several other important changes in the post-shock region: Photons and electron–positron pairs become the dominant source of pressure in the immediate post-shock region, deleptonisation behind the shock occurs more gradually, and the electron neutrino and anti-neutrino luminosities become similar. Most notably, a region of net neutrino heating (gain region) emerges behind the shock. In the ‘delayed neutrino-driven mechanism’ as conceived by Bethe & Wilson (1985) and Wilson (1985), the neutrino heating eventually leads to a sufficient increase of the post-shock pressure to ‘revive’ the shock and make it re-expand, although the post-shock velocity initially remains negative. Since shock expansion increases the mass of the dissociated material exposed to strong neutrino heating, this is thought to be a self-sustaining runaway process that eventually pumps sufficient energy into the post-shock region to allow for the development of positive post-shock velocities and, further down the road, the expulsion of the stellar envelope.

Modern simulations of CCSNe that include energy-dependent neutrino transport, state-of-the art microphysics, and (to various degrees) general relativistic effects have demonstrated that the neutrino-driven mechanism is not viable in spherical symmetry (Rampp & Janka 2000, 2002; Liebendörfer et al. 2001, 2004, 2005; Sumiyoshi et al. 2005; Buras et al. 2006a, 2006b; Müller, Janka, & Dimmelmeier 2010; Fischer et al. 2010; Lentz et al. 2012a, 2012b), except for supernova progenitors of the lowest masses (Kitaura, Janka, & Hillebrandt 2006; Janka et al. 2008; Burrows, Dessart, & Livne 2007b; Fischer et al. 2010), which will be discussed in Section 2.

In its modern guise, the paradigm of neutrino-driven explosions therefore relies on the joint action of neutrino heating and various hydrodynamic instabilities to achieve shock revival. As demonstrated by the first generation of multi-dimensional supernova models in the 1990s (Herant et al. 1994; Burrows, Hayes, & Fryxell 1995; Janka & Müller 1995, 1996), the gain region is subject to convective instability due to the negative entropy gradient established by neutrino heating. Convection can be suppressed if the accreted material is quickly advected from the shock to the gain radius (Foglizzo, Scheck, & Janka 2006). Under these conditions, the standing accretion shock instability (SASI; Blondin, Mezzacappa, & DeMarino 2003; Blondin & Mezzacappa 2006; Foglizzo et al. 2007; Laming 2007; Yamasaki & Yamada 2007; Fernández & Thompson 2009a, 2009b) can still grow, which is mediated by an advective-acoustic cycle (Foglizzo 2002; Foglizzo et al. 2007; Guilet & Foglizzo 2012) and manifests itself in the form of large-scale sloshing and spiral motions of the shock. The precise mechanism whereby these instabilities aid shock revival requires careful discussion (see Section 3.3), but their net effect can be quantified using the concept of the ‘critical luminosity’ (Burrows & Goshy 1993) for the transition from a steady-state accretion flow to runaway shock expansion: In effect, convection and/or the SASI reduce the critical luminosity in multi-D by 20. . .30% (Murphy & Burrows 2008; Nordhaus et al. 2010; Hanke et al. 2012; Fernández 2015) compared to the case of spherical symmetry (1D).

1.2 Current questions and structure of this review

We cannot hope to comprehensively review all aspects of the CCSN explosion problem, even if we limit ourselves to the neutrino-driven paradigm. Instead we shall focus on the following topics that immediately connect to the above overview of the neutrino-driven mechanism:

• • The neutrino-driven explosion mechanism demonstrably works at the low-mass end of supernova progenitors. In Section 2, we shall discuss the specific explosion dynamics in the region around the mass limit for iron-core formation, i.e., for ECSN progenitors and structurally similar iron-core progenitors. We shall also consider the nucleosynthesis in these explosions; since they are robust, occur early after bounce, and can easily be simulated until the explosion energy has saturated, explosions of ECSN and ECSN-like progenitors currently offer the best opportunity to study CCSN nucleosynthesis based on first-principle explosion models.

• • For more massive progenitors, it has yet to be demonstrated that the neutrino-driven mechanism can produce robust explosions in 3D with explosion properties (e.g., explosion energy, nickel mass, remnant mass) that are compatible with observations. In Section 3, we shall review the current status of 3D supernova simulations, highlighting the successes and problems of the current generation of models and detailing the recent progress towards a quantitative understanding of the interplay of neutrino heating and multi-dimensional fluid flow.

• • In the wake of a rapid expansion of the field of CCSN modelling, a wide variety of methods have been employed to investigate the supernova problem with a continuum from a rigorous first-principle approach to parameterised models of limited applicability that are only suitable for attacking well-circumscribed problems. In Section 4, we present an overview of the different numerical approaches to simulations of neutrino-driven explosions and provide some guidance for assessing and comparing simulation results.

• • The problem of shock revival by the neutrino-driven mechanism has not been conclusively solved. In Section 5, we shall review one of the promising ideas that could help explain supernova explosions over a wide range of progenitors, viz. the suggestion that shock revival may be facilitated by strong seed perturbations from prior convective shell burning in the infalling O or Si shells (Arnett & Meakin 2011; Couch & Ott 2013; Müller & Janka 2015; Couch et al. 2015; Müller et al. 2016a); and we shall also discuss some other perspectives opened up by current and future 3D simulations of late burning stages in supernova progenitors.

Potential observational probes for multi-dimensional fluid flow in the supernova core during the first $\mathord {\sim } 1 \, {\rm s}$ exist in the form of the neutrino and gravitational wave signals, but we shall not touch these in any depth and instead point the reader to topical reviews (Ott 2009 and Kotake 2013 for gravitational wave emission; Mirizzi et al. 2016 for the neutrino signal) as well as some of the major publications of recent years (gravitational waves: Müller, Janka, & Marek 2013; Yakunin et al. 2015; Nakamura et al. 2016; neutrinos: Tamborra et al. 2013, 2014a; Müller & Janka 2014). Neither do we address alternative explosion scenarios here and refer the reader to Janka (2012) for a broader discussion that covers, e.g., the magnetorotational mechanism as the most likely explanation for hypernovae with explosion energies of up to $\mathord {\sim } 10^{52} \, {\rm erg}$ .

2 THE LOW-MASS END ELECTRON-CAPTURE SUPERNOVAE AND THEIR COUSINS

Stars with ZAMS masses in the range $\mathord {\sim } 8 \ldots 10\,\text{M}_\odot$ exhibit structural peculiarities during their evolution that considerably affect the supernova explosion dynamics if they undergo core-collapse. The classical path towards ECSNe (Nomoto 1984, 1987), where electron captures on ²⁴Mg and ²⁰Ne in a degenerate O–Ne–Mg core of $\mathord {{\sim }1.37}\,\text{M}_\odot$ drive the core towards collapse, best exemplifies these peculiarities: Only a small C/O layer is present on top of the core, and the He layer has been effectively whittled down by dredge-up. The consequence is an extremely steep density gradient between the core and the high-entropy hydrogen envelope (Figure 1). Whilst this particular scenario is beset with many uncertainties (Siess 2007; Poelarends et al. 2008; Jones et al. 2013, 2014, 2016a; Doherty et al. 2015; Schwab et al. 2015; Woosley & Heger 2015b), recent studies of stellar evolution in the mass range around $9\,\text{M}_\odot$ have demonstrated that there is a variety of paths towards core-collapse that result in a similar progenitor structure (Jones et al. 2013; Woosley & Heger 2015b), though there is some variation, e.g., in the mass of the remaining He shell due to a different history of dredge-up events. From the perspective of supernova explosion dynamics, the crucial features in the mass range around $9\,\text{M}_\odot$ are the small mass of the remaining C/O shell and the rapid drop of the density outside the core; both are shared by ECSN progenitors and the lowest iron-core progenitors. This is illustrated in Figure 1 (see also Figure 7 in Jones et al. 2013 and Figure 4 in Woosley & Heger 2015b).

Figure 1. Density profiles of several low-mass supernova progenitors illustrating the conditions for ECSN-like explosions. Profiles are shown for the $8.8\,\text{M}_\odot$ ECSN-progenitor of Nomoto (1984, 1987) (N8.8, black), the 8.8 $\text{M}_\odot$ ‘failed massive star’ of Jones et al. (2013) (J8.8, purple), low-mass iron-core progenitors (A. Heger, private communication) of $9.6\,\text{M}_\odot$ (z9.6, with Z=0, red) and $8.1\,\text{M}_\odot$ (u8.1, with Z = 10⁻⁴, blue), and iron progenitors with $10.09\,\text{M}_\odot$ and $15\,\text{M}_\odot$ (s10.09 and s15, from Müller et al. 2016b, yellow and cyan), and $11.2\,\text{M}_\odot$ (s11.2 from Woosley et al. 2002, green). The thick dashed vertical line roughly denotes the location of the shell that reaches the shock 0.5 s after the onset of collapse. Slanted dashed lines roughly demarcate the regime where the accretion rate onto the shock reaches $0.05\,\text{M}_\odot \, {\rm s}^{-1}$ (thick dashed line), $5 \times 10^{-3}\,\text{M}_\odot \, {\rm s}^{-1}$ (thin), and $5 \times 10^{-4}\,\text{M}_\odot \, {\rm s}^{-1}$ (thin) (see Section 2.1.2 for details and underlying assumptions). ECSN-like explosion dynamics is expected if the density profile intersects the grey region.

2.1 Explosion dynamics in ECSN-like progenitors

2.1.1 Classical electron-capture supernova models

The steep density gradient outside the core in ECSN-like progenitors is immediately relevant for the dynamics of the ensuing supernova because it implies a rapid decline of the mass accretion rate $\dot{M}$ as the edge of the core reaches the stalled accretion shock. A rapid drop in $\dot{M}$ implies a decreasing ram pressure ahead of the shock and a continuously increasing shock radius (though the shock remains a stationary accretion shock for at least $\mathord {\sim } 50\, {\rm ms}$ after bounce and longer for some ECSN-like progenitor models). Under these conditions, neutrino heating can easily pump sufficient energy into the gain region to make the accreted material unbound and power runaway shock expansion. As a result, the neutrino-driven mechanism works for ECSN-like progenitors even under the assumption of spherical symmetry. Using modern multi-group neutrino transport, this was demonstrated by Kitaura et al. (2006) for the progenitor of Nomoto (1984, 1987) and confirmed in subsequent simulations by different groups (Janka et al. 2008; Burrows et al. 2007b; Fischer et al. 2010). The explosions are characterised by a small explosion energy of $\mathord {\sim } 10^{50} \, {\rm erg}$ (Kitaura et al. 2006; Janka et al. 2008) and a small nickel mass of a few $ 10^{-3}\,\text{M}_\odot$ (Wanajo et al. 2009).

Even though multi-dimensional effects are not crucial for shock revival in these models, they are not completely negligible. Higher entropies at the bottom of the gain layer lead to convective overturn driven by Rayleigh–Taylor instability shortly after the explosion is initiated (Wanajo, Janka, & Müller 2011). Simulations in axisymmetry (2D) showed that this leads to a modest increase of the explosion energy in Janka et al. (2008); an effect which is somewhat larger in more recent models (von Groote et al., in preparation). The effect of Rayleigh–Taylor overturn on the ejecta composition is, however, much more prominent (see Section 2.2).

2.1.2 Conditions for ECSN-like explosion dynamics

Not all of the newly available supernova progenitor models at the low-mass end (Jones et al. 2013, 2014; Woosley & Heger 2015b) exhibit a similarly extreme density profile as the model of Nomoto (1984, 1987); in some of them, the density gradient is considerably more shallow (Figure 1). This prompts the questions: How steep a density gradient is required outside the core to obtain an explosion that is triggered by a rapid drop of the accretion rate and works with no or little help from multi-D effects? In reality, there will obviously be a continuum between ECSN-like events and neutrino-driven explosions of more massive stars, in which multi-D effects are crucial for achieving shock revival. Nonetheless, a rough distinction between the two different regimes is still useful, and can be based on the concept of the critical neutrino luminosity of Burrows & Goshy (1993).

Burrows & Goshy (1993) showed that stationary accretion flow onto a proto-neutron star in spherical symmetry is no longer possible if the neutrino luminosity L _ν (which determines the amount of heating) exceeds a critical value $L_{\rm crit}(\dot{M})$ that is well approximated by a power law in $\dot{M}$ with a small exponent, or, equivalently, if $\dot{M}$ drops below a threshold value for a given luminosity. This concept has recently been generalised (Janka 2012; Müller & Janka 2015; Summa et al. 2016; Janka, Melson, & Summa 2016) to a critical relation for the (electron-flavour) neutrino luminosity L _ν and neutrino mean energy E _ν as a function of mass accretion rate $\dot{M}$ and proto-neutron star mass M as well as additional correction factors, e.g., for shock expansion due to non-radial instabilities.

For low-mass progenitors with tenuous shells outside the core, M, L _ν, and E _ν do not depend dramatically on the stellar structure outside the core during the early post-bounce: The proto-neutron star mass is inevitably $M\approx 1.4\,\text{M}_\odot$ , and since the neutrino emission is dominated by the diffusive neutrino flux from the core, the neutrino emission properties are bound to be similar to the progenitor of Nomoto (1984), i.e., one has L _ν ~ 5 × 10⁵² erg s⁻¹ and E _ν ≈ 11 MeV (Hüdepohl et al. 2009), with a steady decrease of the luminosity towards later times. Using calibrated relations for the ‘heating functional’ ² L _ν E ² _ν (Janka et al. 2016), this translates into a critical mass accretion rate of $\dot{M}_{\rm crit} \approx 0.07\,\text{M}_\odot \, {\rm s}^{-1}$ for ECSN-like progenitors.

To obtain similarly rapid shock expansion as for the $8.8\,\text{M}_\odot$ model of Nomoto (1984), $\dot{M}$ must rapidly plummet well below this value. This can be translated into a condition for the density profile outside the core using analytic expressions for the infall time t _infall and accretion rate $\dot{M}$ for mass shell m, which are roughly given by (Woosley & Heger 2012, 2015a; Müller et al. 2016b),

(1) $$\begin{equation} t_{\rm infall}=\sqrt{\frac{\pi }{4 G \bar{\rho }}} =\sqrt{\frac{\pi ^2 r^3}{3 G m}}, \end{equation}$$

and

(2) $$\begin{equation} \dot{M}=\frac{2 m }{t_{\rm infall}} \frac{\rho }{\bar{\rho } - \rho }, \end{equation}$$

where $\bar{\rho }$ is the average density inside the mass shell. For progenitors with little mass outside the core, we have

(3) $$\begin{equation} \dot{M}\approx \frac{2 m }{t_{\rm infall}} \frac{\rho }{\bar{\rho }} =\frac{8\rho }{3}\sqrt{3 G m r^3}. \end{equation}$$

Using $m=1.4\,\text{M}_\odot$ and assuming that $\dot{M}$ needs to drop at least to $M_{\rm crit}=0.05\,\text{M}_\odot \, {\rm s}^{-1}$ within 0.5 s after the onset of collapse to obtain ECSN-like explosion dynamics, one finds that the density needs to drop to

(4) $$\begin{equation} \rho \lesssim \frac{1}{8}\sqrt{\frac{3}{G m}}\dot{M}_{\rm crit} r^{-3/2} \end{equation}$$

for a radius r < 2 230 km.

Figure 1 illustrates that the density gradient at the edge of the core can be far less extreme than in the model of Nomoto (1984) to fulfil this criterion. ECSN-like explosion dynamics is expected alike for the modern $8.8\,\text{M}_\odot$ ECSN progenitor of Jones et al. (2013) and low-mass iron cores (A. Heger, private communication) of $8.1\,\text{M}_\odot$ (with metallicity Z = 10⁻⁴) and $9.6\,\text{M}_\odot$ (Z = 0), though the low-mass iron-core progenitors are a somewhat marginal case.

2.1.3 Low-mass iron-core progenitors

Simulations of these two low-mass iron progenitors with $8.1\,\text{M}_\odot$ (Müller, Janka, & Heger 2012b) and $9.6\,\text{M}_\odot$ (Janka et al. 2012; Müller et al. 2013 in 2D; Melson, Janka, & Marek 2015a in 3D) nonetheless demonstrated that the structure of these stars is sufficiently extreme to produce explosions reminiscent of ECSN models: Shock revival sets in early around 100 ms after bounce, aided by the drop of the accretion rate associated with the infall of the thin O and C/O shells, and the explosion energy remains small (5 × 10⁴⁹. . .10⁵⁰ erg).

As shown by Melson et al. (2015a), there are important differences to ECSNe, however: Whilst shock revival also occurs in spherical symmetry, multi-dimensional effects significantly alter the explosion dynamics. In 1D, the shock propagates very slowly through the C/O shell after shock revival, and only accelerates significantly after reaching the He shell. Without the additional boost by convective overturn, the explosion energy is lower by a factor of $\mathord {\sim } 5$ compared to the multi-D case. Different from ECSNe, somewhat slower shock expansion provides time for the small-scale convective plumes to merge into large structures as shown for the $9.6\,\text{M}_\odot$ model of Janka et al. (2012) in Figure 2.

Figure 2. Entropy s (left half of plot) and electron fraction $Y_\text{e}$ (right half) in the $9.6\,\text{M}_\odot$ explosion model of Janka et al. (2012) and Müller et al. (2013) 280 ms after bounce. Large convective plumes push neutron-rich material from close to the gain region out at high velocities.

Both for the $8.8\,\text{M}_\odot$ model of Wanajo et al. (2011) and the low-mass iron-core explosion models, the dynamics of the Rayleigh–Taylor plumes developing after shock revival is nonetheless quite similar. The entropy of the rising plumes is roughly $\mathord {\sim } 15 \ldots 20 k_{\rm b}/{\rm nucleon}$ compared to $\mathord {\sim }10 k_{\rm b}/{\rm nucleon}$ in the ambient medium. For such an entropy contrast, balance between buoyancy and drag forces applies a limiting velocity of the order of the speed of sound. This limit appears to be reached relatively quickly in the simulations. Apart from the very early growth phase, the plume velocities should therefore not depend strongly on the initial seed perturbations; they are rather set by bulk parameters of the system, namely the post-shock entropy at a few hundred kilometres and the entropy close to the gain radius, which together determine the entropy contrast of the plumes. This will become relevant later in our discussion of the nucleosynthesis of ECSN-like explosions.

2.2 Nucleosynthesis

2.2.1 1D electron-capture supernovae models—early ejecta

Nucleosynthesis calculations based on modern, spherically symmetric ECSN models were first performed by Hoffman, Müller, & Janka (2008) and Wanajo et al. (2009). The results of these calculations appeared to point to a severe conflict with observational constraints, showing a strong overproduction of N = 50 nuclei, in particular ⁹⁰Zr, due to the ejection of slightly neutron-rich material (electron fraction $Y_\text{e}\gtrsim 0.46$ ) with relatively low entropy (s ≈ 18k_b /nucleon) immediately after shock revival. Hoffman et al. (2008) inferred that such nucleosynthesis yields would only be compatible with chemogalactic evolution if ECSNe were rare events occurring at a rate no larger than once per 3 000 yr.

The low $Y_\text{e}$ -values in the early ejecta stem from the ejection of matter at relatively high velocities in the wake of the fast-expanding shock. In slow outflows, neutrino absorption on neutrons and protons drives $Y_\text{e}$ to an equilibrium value that is set by the electron neutrino and anti-neutrino luminosities $L_{\nu _\text{e}}$ and $L_{\bar{\nu }_\text{e}}$ , the ‘effective’ mean energies ³ $\varepsilon _{\nu _\text{e}}$ and $\varepsilon _{{\bar{\nu }}_\text{e}}$ , and the proton–neutron mass difference Δ = 1.293 MeV as follows (Qian & Woosley 1996):

(5) $$\begin{equation} Y_\text{e} \approx \left[1+\frac{L_{\bar{\nu }_\text{e}} (\varepsilon _{\bar{\nu }_\text{e}}-2 \Delta )}{L_{{\nu }_\text{e}} (\varepsilon _{{\nu }_\text{e}}+2 \Delta )}\right]^{-1}. \end{equation}$$

For the relatively similar electron neutrino and anti-neutrino luminosities and a small difference in the mean energies of 2. . .3 MeV in modern simulations, one typically finds an asymptotic value of $Y_\text{e}>0.5$ , i.e., proton-rich conditions. To obtain low $Y_\text{e}<0.5$ in the ejecta, neutrino absorption reactions need to freeze out at a high density (small radius) when the equilibrium between the reactions $n(\nu _\text{e},\text{e}^-)p$ and $p(\nu _\text{e},\text{e}^+)n$ is still skewed towards low $Y_\text{e}$ due to electron captures $p (\text{e}^-,\nu _\text{e}) n$ on protons. Neglecting the difference between arithmetic, quadratic, and cubic neutrino mean energies and assuming a roughly equal contribution of $n(\nu _\text{e},\text{e}^-)p$ and $p(\bar{\nu }_\text{e},\text{e}^+)n$ to the neutrino heating, one can estimate that freeze-out roughly occurs when (cp. Equation (81) in Qian & Woosley 1996),

(6) $$\begin{equation} \frac{v_\text{r}}{r} \approx \frac{2 m_\text{N} \dot{q}_\nu }{E_{\nu _\text{e}}+E_{\bar{\nu }_\text{e}}}, \end{equation}$$

where $m_\text{N}$ is the nucleon mass, $\dot{q}_\nu$ is the mass-specific neutrino heating rate, r is the radius, and $v_\text{r}$ is the radial velocity. Since $\dot{q}_\nu \propto r^{-2}$ , freeze-out will occur at smaller r, higher density, and smaller $Y_\text{e}$ for higher ejection velocity.

2.2.2 Multi-D effects and the composition of the early ejecta

Since high ejection velocities translate into lower $Y_\text{e}$ , the Rayleigh–Taylor plumes in 2D simulations of ECSNe (Figure 2 in Wanajo et al. 2011) and explosions of low-mass iron cores (Figure 2) contain material with even lower $Y_\text{e}$ than found in 1D ECSN models. Values of $Y_\text{e}$ as low as 0.404 are found in Wanajo et al. (2011).

Surprisingly, Wanajo et al. (2011) found that the neutron-rich plumes did not aggravate the problematic overproduction of N = 50 nuclei in their 2D ECSN model. This is due to the fact that the entropy in the neutron-rich lumps is actually smaller than in 1D ⁴ (but higher than in the ambient medium), which changes the character of the nucleosynthesis by reducing the α-fraction at freeze-out from nuclear statistical equilibrium (NSE). The result is an interesting production of trans-iron elements between Zn and Zr for the progenitor of Nomoto (1984, 1987); the production factors are consistent with current rate estimates for ECSNe of about 4% of all supernovae (Poelarends et al. 2008). Subsequent studies showed that neutron-rich lumps in the early ejecta of ECSNe could contribute a sizeable fraction to the live ⁶⁰Fe in the Galaxy (Wanajo, Janka, & Müller 2013b), and might be production sites for some other rare isotopes of obscure origin, such as ⁴⁸Ca (Wanajo, Janka, & Müller 2013a). Due to the similar explosion dynamics, low-mass iron-core progenitors exhibit rather similar nucleosynthesis (Wanajo et al., in preparation; Harris et al., in preparation). The results of these nucleosynthesis calculations tallies with the observed abundance trends in metal-poor stars that suggest a separate origin of elements like Sr, Y, and Zr from the heavy r-process elements (light element primary process; Travaglio et al. 2004; Wanajo & Ishimaru 2006; Qian & Wasserburg 2008; Arcones & Montes 2011; Hansen et al. 2012; Ting et al. 2012).

Since $Y_\text{e}$ in the early ejecta of ECSNe and ECSN-like explosion is sensitive to the neutrino luminosities and mean energies and to the ejection velocity of the convective plumes (which may be different in 3D compared to 2D, or exhibit stochastic variations), Wanajo et al. (2011) also explored the effect of potential uncertainties in the minimum $Y_\text{e}$ in the ejecta on the nucleosynthesis. They found that a somewhat lower $Y_\text{e}$ of ~ 0.3 in the plumes might make ECSNe a site for a ‘weak r-process’ that could explain the enhanced abundances of lighter r-process elements up to Ag and Pd in some metal-poor halo stars (Wanajo & Ishimaru 2006; Honda et al. 2006).

Whether the neutron-rich conditions required for a weak r-process can be achieved in ECSNe or low-mass iron-core supernovae remains to be determined. Figure 3 provides a tentative glimpse on the effects of stochasticity and dimensionality on the $Y_\text{e}$ in neutron-rich plumes based on several 2D and 3D explosion models of a $9.6\,\text{M}_\odot$ low-mass iron-core progenitor (A. Heger, private communication) conducted using the FMT transport scheme of Müller & Janka (2015). ⁵ Stochastic variations in 2D models due to different (random) initial perturbations shift the minimum $Y_\text{e}$ in the ejecta at most by 0.02. This is due to the fact that the Rayleigh–Taylor plumes rapidly transition from the initial growth phase to a stage where buoyancy and drag balance each other and determine the velocity (Alon et al. 1995). 3D effects do not change the distribution of $Y_\text{e}$ tremendously either, at best they tend to shift it to slightly higher values compared to 2D, which is consistent with a somewhat stronger braking of expanding bubbles in 3D as a result of the forward turbulent cascade (Melson et al. 2015a). It thus appears unlikely that the dynamics of convective overturn is a major source of uncertainty for the nucleosynthesis in ECSN-like explosions, though confirmation with better neutrino transport is still needed.

Figure 3. Binned distribution of the electron fraction $Y_\text{e}$ in the early ejecta for different explosion models of a $9.6\,\text{M}_\odot$ star 270 ms after bounce. The plots show the relative contribution ΔM _ej/M _ej to the total mass of (shocked) ejecta in bins with $\Delta Y_\text{e}=0.01$ . The upper panel shows the $Y_\text{e}$ -distribution for the 2D model of Janka et al. (2012) computed using the vertex-coconut code (Müller et al. 2010). The bottom panel illustrates the effect of stochastic variations and dimensionality using several 2D models (thin lines) and a 3D model computed with the coconut-fmt code of Müller & Janka (2015) (thick lines). Note that the dispersion in $Y_\text{e}$ in the early ejecta is similar for both codes, though the average $Y_\text{e}$ in the early ejecta is spuriously low when less accurate neutrino transport is used (fmt instead of Vertex). The bottom panel is therefore only intended to show differential effects between different models, and is not a prediction of the absolute value of $Y_\text{e}$ . It suggests that (i) stochastic variations do not strongly affect the distribution of $Y_\text{e}$ in the ejecta, and that (ii) the resulting distribution of $Y_\text{e}$ in 2D and 3D is relatively similar.

If these events are indeed sites of a weak r-process, the missing ingredient is likely to be found elsewhere. Improvements in the neutrino opacities, such as the proper inclusion of nucleon potentials in the charged-current interaction rates (Martínez-Pinedo et al. 2012; Roberts, Reddy, & Shen 2012), or flavour oscillations involving sterile neutrinos (Wu et al. 2014) could lower $Y_\text{e}$ somewhat. Wu et al. (2014) found a significant reduction of $Y_\text{e}$ by up to 0.15 in some of the ejecta, but these results may depend sensitively on the assumption that collective flavour oscillations are still suppressed during the phase in question. Moreover, Wu et al. (2014) pointed out that a reduction of $Y_\text{e}$ with the help of active-sterile flavour conversion might require delicate fine-tuning to avoid shutting off neutrino heating before the onset of the explosion due to the disappearance of $\nu _\text{e}$ ’s (which could be fatal to the explosion mechanism).

Moreover, whether ECSNe necessarily need to co-produce Ag and Pd with Sr, Y, and Zr is by no means clear. Whilst observed abundance trends may suggest such a co-production, the abundance patterns of elements between Sr and Ag in metal-poor stars appear less robust (Hansen, Montes, & Arcones 2014); and the failure of unaltered models to produce Ag and Pd may not be indicative of a severe tension with observations.

2.2.3 Other nucleosynthesis scenarios for electron-capture supernovae

There are at least two other potentially interesting sites for nucleosynthesis in ECSN-like supernovae. For ‘classical’ ECSN-progenitors with more extreme density profiles, it has been proposed that the rapid acceleration of the shock in the steep density gradient outside the core can lead to sufficiently high post-shock entropies (s ~ 100 k_b /nucleon) and short expansion time-scales (τ_exp ~ 10⁻⁴ s) to allow r-process nucleosynthesis in the thin shells outside the core (Ning, Qian, & Meyer 2007). This has not been borne out by numerical simulations, however (Janka et al. 2008; Hoffman et al. 2008). When the requisite high entropy is reached, the post-shock temperature has already dropped far too low to dissociate nuclei, and the expansion timescale does not become sufficiently short for the scenario of Ning et al. (2007) to work. The proposed r-process in the rapidly expanding shocked shells would require significantly different explosion dynamics, e.g., a much higher explosion energy.

The neutrino-driven wind that is launched after accretion onto the proto-neutron star has been completely subsided has long been discussed as a potential site of r-process nucleosynthesis in supernovae (Woosley et al. 1994; Takahashi, Witti, & Janka 1994; Qian & Woosley 1996; Cardall & Fuller 1997; Thompson, Burrows, & Meyer 2001; Arcones, Janka, & Scheck 2007; Arcones & Thielemann 2013). ECSN-like explosions are in many respects the least favourable site for an r-process in the neutrino-driven wind since they produce low-mass neutron stars, which implies low wind entropies and long expansion timescales (Qian & Woosley 1996), i.e., conditions that are detrimental to r-process nucleosynthesis. However, ECSNe are unique in as much as the neutrino-driven wind can be calculated self-consistently with Boltzmann neutrino transport (Hüdepohl et al. 2009; Fischer et al. 2010) without the need to trigger an explosion artificially. These simulations revealed a neutrino-driven wind that is not only of moderate entropy (s ≲ 140k_b /nucleon even at late times), but also becomes increasingly proton-rich with time, in which case the νp-process (Fröhlich et al. 2006) could potentially operate. The most rigorous nucleosynthesis calculations for the neutrino-driven wind in ECSNe so far (Pllumbi et al. 2015) are based on simulations that properly account for nucleon interaction potentials in the neutrino opacities (Martínez-Pinedo et al. 2012; Roberts et al. 2012) and have also explored the effects of collective flavour oscillations, active-sterile flavour conversion. Pllumbi et al. (2015) suggest that wind nucleosynthesis in ECSNe is rather mundane: Neither does the νp-process operate nor can neutron-rich conditions be restored to obtain conditions even for a weak r-process. Instead, they find that wind nucleosynthesis mainly produces nuclei between Sc and Zn, but the production factors are low, implying that the role of neutrino-driven winds in ECSNe is negligible for this mass range for the purpose of chemogalactic evolution.

2.3 Electron-capture supernovae—transients and remnants

Although the explosion mechanism of ECSNe is in many respects best understood amongst all CCSN types from the viewpoint of explosion mechanism, unambiguously identifying transients as ECSNe has proved more difficult. It has long been proposed that SN 1054 was an ECSN (Nomoto et al. 1982) based on the properties of its remnant, the Crab nebula: The total mass of ejecta in the nebula is small ( ${\lesssim }5\,\text{M}_\odot$ ; Davidson & Fesen 1985; MacAlpine & Uomoto 1991; Fesen, Shull, & Hurford 1997), as is the oxygen abundance (Davidson et al. 1982; Henry & MacAlpine 1982; Henry 1986), which is in line with the thin O-rich shells in ECSN progenitors. Moreover, the kinetic energy of the ejecta is only about ≲ 10⁵⁰ erg (Fesen et al. 1997; Hester 2008) as expected for an ECSN-like event. Whether the Crab originates from a classical ECSN or from something slightly different like a ‘failed massive star’ of Jones et al. (2013) continues to be debated; MacAlpine & Satterfield (2008) have argued, for example, against the former interpretation based on a high abundance ratio of C vs. N and the detection of some ashes of oxygen burning (S, Ar) in the nebula.

It has been recognised in recent years that the (reconstructed) light curve of SN 1054—a type IIP supernova with a relatively bright plateau—is also compatible with the low explosion energy of ≲ 10⁵⁰ erg predicted by recent numerical simulations. Smith (2013) interpreted the bright plateau, which made SN 1054 visible by daytime for $\mathord {\sim } 3 \ {\rm weeks}$ , as the result of interaction with circumstellar medium (CSM). The scenario of Smith (2013) requires significant mass loss ( $0.1\,\text{M}_\odot$ for about 30 yr) shortly before the supernova, which may be difficult to achieve, although some channels towards ECSN-like explosions could involve dramatic mass loss events (Woosley & Heger 2015b). Subsequent numerical calculations of ECSN light curves (Tominaga, Blinnikov, & Nomoto 2013; Moriya et al. 2014) demonstrated, however, that less extreme assumptions for the mass loss are required to explain the optical signal of SN 1054; indeed a very extended hydrogen envelope may be sufficient to explain the bright plateau, and CSM interaction with the progenitor wind may only be required to prevent the SN from fading too rapidly.

Several other transients have also been interpreted as ECSNe, e.g., faint type IIP supernovae such as SN 2008S (Botticella et al. 2009). Smith (2013) posits that ECSNe are observed type IIn-P supernovae with circumstellar interaction like SN 1994W with a bright plateau and a relatively sharp drop to a faint nickel-powered tail, but again the required amount of CSM is not easy to explain. All of these candidate events share low kinetic energies and small nickel masses as a common feature and are thus prima facie compatible with ECSN-like explosion dynamics. Variations in the envelope structure of ECSN-progenitors (e.g., envelope stripping in binaries) may account for the very different optical signatures (Moriya et al. 2014).

The peculiar nucleosynthesis in ECSNe-like explosions may also leave observable fingerprints in the electromagnetic signatures. The slightly neutron-rich character of the early ejecta results in a strongly supersolar abundance ratio of Ni to Fe after β-decays are completed (Wanajo et al. 2011). Such high Ni/Fe ratios are seen in the nebular spectra of some supernovae (Jerkstrand et al. 2015a, 2015b). ECSNe can only explain some of these events; however, many of them exhibit explosion energies and Nickel masses that are incompatible with an ECSN.

3 3D SUPERNOVA MODELS OF MASSIVE PROGENITORS

In more massive progenitors with extended Si and O shells, the mass accretion rate onto the shock does not drop as rapidly as in ECSN-like explosions. Typically, one finds a relatively stable accretion rate of a few $0.1\,\text{M}_\odot \, {\rm s}^{-1}$ during the infall of the O shell, which implies a high ram pressure ahead of the shock. Under these conditions, it is no longer trivial to demonstrate that neutrino heating can pump a sufficient amount of energy into the post-shock region to power runaway shock expansion. 1D simulations of the post-bounce phase using Boltzmann solvers for the neutrino transport convincingly demonstrated that neutrino-driven explosions cannot be obtained under such conditions in spherical symmetry (Liebendörfer et al. 2001; Rampp & Janka 2000; Burrows et al. 2000a). Much of the work of recent years has therefore focussed on better understanding and accurately modelling how multi-dimensional effects in supernovae facilitate neutrino-driven explosions—an undertaking first begun in the 1990s with axisymmetric (2D) simulations employing various approximations for neutrino heating and cooling (Herant, Benz, & Colgate 1992; Yamada, Shimizu, & Sato 1993; Herant et al. 1994; Burrows et al. 1995; Janka & Müller 1995, 1996). 2D simulations have by now matured to the point that multi-group neutrino transport and the neutrino-matter interactions can be modelled with the same rigour as in spherical symmetry (Livne et al. 2004; Buras et al. 2006a; Müller et al. 2010; Bruenn et al. 2013; Just, Obergaulinger, & Janka 2015; Skinner, Burrows, & Dolence 2016), or with, still with acceptable accuracy for many purposes (see Section 4 for a more careful discussion), by using some approximations either in the transport treatment or the neutrino microphysics (Suwa et al. 2010; Müller & Janka 2015; Pan et al. 2016; O’Connor & Couch 2015; Roberts et al. 2016).

3.1 Prelude—first-principle 2D models

The current generation of 2D supernova simulations with multi-group neutrino transport has gone a long way towards demonstrating that neutrino heating can bring about explosion in conjunction with convection or the SASI. Thanks to steadily growing computational resources, the range of successful neutrino-driven explosion models has grown from about a handful in mid-2012 (Buras et al. 2006b; Marek & Janka 2009; Suwa et al. 2010; Müller, Janka, & Marek 2012a) to a huge sample of explosion models with ZAMS masses between 10 and $ 75\,\text{M}_\odot$ , different metallicities, and different choices for the supranuclear equation of state (Müller et al. 2012b; Janka et al. 2012; Suwa et al. 2013; Bruenn et al. 2013; Obergaulinger, Janka, & Aloy 2014; Nakamura et al. 2015; Müller 2015; Bruenn et al. 2016; O’Connor & Couch 2015; Summa et al. 2016; Pan et al. 2016).

Many of the findings from these simulations remain important and valid after the advent of 3D modelling: The 2D models have established, amongst other things, the existence of distinct SASI- and convection-dominated regimes in the accretion phase, both of which can lead to successful explosion (Müller et al. 2012b) in agreement with tunable, parameterised models (Scheck et al. 2008; Fernández et al. 2014). They have shown that ‘softer’ nuclear equations of state that result in more compact neutron stars are generally favourable for shock revival (Janka 2012; Suwa et al. 2013; Couch 2013a). The inclusion of general relativistic effects, whether by means of the conformally flat approximation (CFC) or, less rigorously, an effective pseudo-relativistic potential for Newtonian hydrodynamics, was found to have a similarly beneficial effect (CFC: Müller et al. 2012a; pseudo-Newtonian: O’Connor & Couch 2015). Moreover, there are signs that the 2D models of some groups converge with each other; simulations of four different stellar models ( $12,15,20,25\,\text{M}_\odot$ ) of Woosley & Heger (2007) by Summa et al. (2016) and O’Connor & Couch (2015) have yielded quantitatively similar results.

Despite these successes, 2D models have, by and large, struggled to reproduce the typical explosion properties of supernovae. They are often characterised by a slow and unsteady growth of the explosion energy after shock revival. Usually the growth of the explosion energy cannot be followed beyond 2. . .4 × 10⁵⁰ erg after simulating up to $\mathord {\sim } 1 \, {\rm s}$ of physical time (Janka et al. 2012; Nakamura et al. 2015; O’Connor & Couch 2015), i.e., below typical observed values of 5. . .9 × 10⁵⁰ erg (Kasen & Woosley 2009; Pejcha & Prieto 2015). Only the models of Bruenn et al. (2016) reach significantly higher explosion energies. Whilst the explosion energy often has not levelled out yet at the end of the simulations and may still grow significantly for several seconds (Müller 2015), its continuing growth comes at the expense of long-lasting accretion onto the proto-neutron star. This may result in inordinately high remnant masses. Thus, whilst 2D models appeared to have solved the problem of shock revival, they faced an energy problem instead.

3.2 Status of 3D core-collapse supernova models

Before 3D modelling began in earnest (leaving aside tentative sallies into 3D by Fryer & Warren 2002), it was hoped that 3D effects might facilitate shock revival even at earlier times than in 2D, and that this might then also provide a solution to the energy problem, since more energy can be pumped into the neutrino-heated ejecta at early times when the mass in the gain region is larger. These hopes were already disappointed once several groups investigated the role of 3D effects in the explosion mechanism using a simple ‘light-bulb’ approach, where the neutrino luminosity and mean energy during the accretion phase are prescribed and very simple approximations for the neutrino heating and cooling terms are employed. Although Nordhaus et al. (2010) initially claimed a significant reduction of the critical neutrino luminosity for shock revival in 3D compared to 2D based on such an approach, these results were affected by the gravity treatment (Burrows, Dolence, & Murphy 2012) and have not been confirmed by subsequent studies. Similar parameterised simulations have shown that the critical luminosity in 3D is roughly equal to 2D (Hanke et al. 2012; Couch 2013b; Burrows et al. 2012; Dolence et al. 2013) and about 20% lower than in 1D, though the results differ about the hierarchy between 2D and 3D.

Subsequent supernova models based on multi-group neutrino transport yielded even more unambiguous results: Shock revival in 3D was either not achieved for progenitors that explode in 2D (Hanke et al. 2013; Tamborra et al. 2014b), or was delayed significantly (Takiwaki, Kotake, & Suwa 2014; Melson et al. 2015b; Lentz et al. 2015). These first disappointing results need to be interpreted carefully, however: A detailed analysis of the heating conditions in the non-exploding 3D models of 11.2, 20, and $27\,\text{M}_\odot$ progenitors simulated by the Garching supernova group revealed that these are very close to shock revival (Hanke et al. 2013; Hanke 2014; Melson et al. 2015b). Moreover, the 3D models of the Garching group are characterised by more optimistic heating conditions, larger average shock radii, and higher kinetic energies in non-spherical motions compared to 2D for extended periods of time; the same is true for the delayed (compared to 2D) 3D explosion of Lentz et al. (2015) of a $15\,\text{M}_\odot$ progenitor. It is merely when it comes to sustaining shock expansion that the 3D models prove less resilient than their 2D counterparts, which transition into an explosive runaway more robustly.

The conclusion that 3D models are only slightly less prone to explosion is reinforced by the emergence of the first successful simulations of shock revival in progenitors with $20\,\text{M}_\odot$ (Melson et al. 2015b) and $15\,\text{M}_\odot$ (Lentz et al. 2015) using rigorous multi-group neutrino transport and the best available neutrino interaction rates. There is also a number of 3D explosion models based on more simplified approaches to multi-group neutrino transport (Takiwaki, Kotake, & Suwa 2012; Takiwaki et al. 2014; Müller 2015; Roberts et al. 2016).

3.3 How do multi-D effects facilitate shock revival?

Despite these encouraging developments, several questions now need to be addressed to make further progress: What is the key to robust 3D explosion models across the entire progenitor mass range for which we observe explosions (i.e., at least up to $15 \ldots 18\,\text{M}_\odot$ ; see Smartt et al. 2009 and Smartt 2015)? This question is tightly connected to another, more fundamental one, namely: What are the conditions for an explosive runaway, and how do multi-dimensional effects modify them?

3.3.1 Conditions for runaway shock expansion

Even without the complications of multi-D fluid flow, the physics of shock revival is subtle. In spherical symmetry, one can show that for a given mass accretion rate $\dot{M}$ , there is a maximum (critical) electron-flavour luminosity L _ν at the neutrinosphere above which stationary accretion flow onto the proto-neutron star is no longer possible (Burrows & Goshy 1993; cp. Section 2). This also holds true if the contribution of the accretion luminosity due to cooling outside the neutrinosphere is taken into account (Pejcha & Thompson 2012). The limit for the existence of stationary solutions does not perfectly coincide with the onset of runaway shock expansion, however. Using 1D light-bulb simulations (i.e., neglecting the contribution of the accretion luminosity), Fernández (2012) and Gabay, Balberg, & Keshet (2015) showed that the accretion flow becomes unstable to oscillatory and non-oscillatory instability slightly below the limit of Burrows & Goshy (1993). Moreover, it is unclear whether the negative feedback of shock expansion on the accretion luminosity and hence on the neutrino heating could push models into a limit cycle (cp. Figure 28 of Buras et al. 2006a) even above the threshold for non-stationarity.

Since an a priori prediction of the critical luminosity, $L_\nu (\dot{M})$ is not feasible, heuristic criteria have been developed (Janka & Keil 1998; Janka, Kifonidis, & Rampp 2001; Thompson 2000; Thompson, Quataert, & Burrows 2005; Buras et al. 2006b; Murphy & Burrows 2008; Pejcha & Thompson 2012; Fernández 2012; Gabay et al. 2015; Murphy & Dolence 2015) to gauge the proximity of numerical supernova models to an explosive runaway (rather than for pinpointing the formal onset of the runaway after the fact, which is of less interest). The most commonly used criticality parameters are based on the ratio of two relevant timescales for the gain region (Janka & Keil 1998; Janka et al. 2001; Thompson 2000; Thompson et al. 2005; Buras et al. 2006b; Murphy & Burrows 2008), namely the advection or dwell time τ_adv that accreted material spends in the gain region, and the heating timescale τ_heat over which neutrino energy deposition changes the total or internal energy of the gain region appreciably. If τ_adv > τ_heat, neutrino heating can equalise the net binding energy of the accreted material before it is lost from the gain region, and one expects that the shock must expand significantly due to the concomitant increase in pressure. Since this expansion further increases τ_adv, an explosive runaway is likely to ensue.

The timescale criterion τ_adv/τ_heat > 1 has the virtue of being easy to evaluate since the two timescales can be defined in terms of global quantities such as the total energy E _{tot, g} in the gain region, the volume-integrated neutrino heating rate $\dot{Q}_\nu$ , and the mass M _g in the gain region (which can be used to define $\tau _{\rm adv}=M_{\rm gain}/\dot{M}$ under steady-state conditions). The significance of these global quantities for the problem of shock revival is immediately intuitive, though care must be taken to define the heating timescale properly. Thompson (2000), Thompson et al. (2005), Murphy & Burrows (2008), and Pejcha & Thompson (2012) define τ_heat as the timescale for changes in the internal energy E _int in the gain region:

(7) $$\begin{equation} \tau _{\rm heat} =\frac{E_{\rm int}}{\dot{Q}_\nu } \end{equation}$$

based on the premise that shock expansion is regulated by the increase in pressure (and hence in internal energy). This definition yields unsatisfactory results, however. The criticality parameter can be spuriously low at shock revival if this definition is used (τ_adv/τ_heat < 0.4).

By defining τ_heat in terms of the total (internal+kinetic+potential) energy ⁶ of the gain region (Buras et al. 2006b):

(8) $$\begin{equation} \tau _{\rm heat} =\frac{E_{\rm tot,g}}{\dot{Q}_\nu }, \end{equation}$$

the criterion τ_adv/τ_heat > 1 becomes a very accurate predictor for non-oscillatory instability (Fernández 2012; Gabay et al. 2015). This indicates that the relevant energy scale to which the quasi-hydrostatic stratification of the post-shock region is the total energy (or perhaps the total or stagnation enthalpy) of the gain region, and not the internal energy. This is consistent with the observation that runaway shock expansion occurs roughly once the total energy or the Bernoulli integral (Fernández 2012; Burrows et al. 1995) reach positive values somewhere (not everywhere) in the post-shock region, which is essentially what the timescale criterion estimates. What is crucial is that the density and pressure gradients between the gain radius and the shock (and hence the shock position) depends sensitively on the ratio of enthalpy h (or the internal energy) and the gravitational potential, rather than on enthalpy alone. Under the (justified) assumption that quadratic terms in $v_\text{r}^2$ in the momentum and energy equation are sufficiently small to be neglected in the post-shock region, one can show (see Appendix A) that the logarithmic derivative of the density ρ in the gain region is constrained by

(9) $$\begin{equation} \frac{\partial \ln \rho }{\partial \ln r} >-\frac{3GM}{rh}, \end{equation}$$

where M is the proto-neutron star mass. Once h > GM/r or even e _int > GM/r (where e _int is the internal energy per unit mass), significant shock expansion must ensue due to the flattening of pressure and density gradients.

Janka (2012), Müller & Janka (2015), and Summa et al. (2016) have also pointed out that the timescale criterion can be converted into a scaling law for the critical electron-flavour luminosity L _ν and mean energy E _ν in terms of the proto-neutron star mass M, the accretion rate $\dot{M}$ , and the gain radius r _g,

(10) $$\begin{equation} (L_\nu E_\nu ^2)_{\rm crit} \propto (\dot{M} M)^{3/5} r_{\rm g}^{-2/5}. \end{equation}$$

The concept of the critical luminosity, the timescale criterion, and the condition of positive total energy or a positive Bernoulli parameter at the gain radius are thus intimately related and appear virtually interchangeable considering that they remain approximate criteria for runaway shock expansion anyway. This is also true for some other explosion criteria that have been proposed, e.g., the antesonic condition of Pejcha & Thompson (2012), which states that the sound speed c _s must exceed a certain fraction of the escape velocity v _esc for runaway shock expansion somewhere in the accretion flow:

(11) $$\begin{equation} c_{\rm s}^2>3/16 v_{\rm esc}^2. \end{equation}$$

Approximating the equation of state as a radiation-dominated gas with an adiabatic index γ = 4/3 and a pressure of P = ρe _int/3 = ρh/4, one finds that the antesonic condition roughly translates to

(12) $$\begin{equation} \frac{c_{\rm s}^2}{3/16 v_{\rm esc}^2} =\frac{4/3 P/\rho }{3/8 GM/r} =\frac{32 e_{\rm int}}{27 GM/r} =\frac{8h}{9GM/r}>1, \end{equation}$$

i.e., the internal energy and the enthalpy must be close to the gravitational binding energy (even if the precise critical values for e _int and h may shift a bit for a realistic equation of state). ⁷

3.3.2 Impact of multi-D effects on the heating conditions

Why do multi-D effects bring models closer to shock revival, and how is this reflected in the aforementioned explosion criteria? Do these explosion criteria even remain applicable in multi-D in the first place?

The canonical interpretation has long been that the runaway condition τ_adv > τ_heat remains the decisive criterion in multi-D, and that multi-D effects facilitate shock revival mainly by increasing the advection time-scale τ_adv (Buras et al. 2006b; Murphy & Burrows 2008). Especially close to criticality, τ_heat is also shortened due to feedback processes—better heating conditions imply that the net binding energy in the gain region and hence τ_heat must decrease.

Whilst simulations clearly show increased advection timescales in multi-D compared to 1D (Buras et al. 2006b; Murphy & Burrows 2008; Hanke et al. 2012) as a result of larger shock radii, the underlying cause for larger accretion shock radii in multi-D is more difficult to pinpoint. Ever since the first 2D simulations, both the transport of neutrino-heated high-entropy material from the gain radius out to the shock (Herant et al. 1994; Janka & Müller 1996) as well as the ‘turbulent pressure’ of convective bubbles colliding with the shock (Burrows et al. 1995) have been invoked to explain larger shock radii in multi-D. Both effects are plausible since they change the components P (thermal pressure) and ρv⊗v (where v is the velocity) of the momentum stress tensor that must balance the ram pressure upstream of the shock during stationary accretion.

That the turbulent pressure plays an important role follows already from the high turbulent Mach number $\mathord {{\sim }0.5}$ in the post-shock region (Burrows et al. 1995; Müller et al. 2012b) before the onset of shock revival, and has been demonstrated quantitatively by Murphy, Dolence, & Burrows (2013) and Couch & Ott (2015) using spherical Reynolds decomposition to analyse parameterised 2D and 3D simulations. Using a simple estimate for the shock expansion due to turbulent pressure, Müller & Janka (2015) were even able to derive the reduction of the critical heating functional in multi-D compared to 1D in terms of the average squared turbulent Mach number ⟨Ma²⟩ in the gain region,

(13) $$\begin{eqnarray} (L_\nu E_\nu ^2)_{\rm crit,2D} &\approx & (L_\nu E_\nu ^2)_{\rm crit,1D} \left(1+\frac{4\langle {\rm Ma}^2 \rangle }{3}\right)^{-3/5}\\ \nonumber &\propto & (\dot{M} M)^{3/5} r_{\rm g}^{-2/5} \left(1+\frac{4\langle {\rm Ma}^2 \rangle }{3}\right)^{-3/5}, \end{eqnarray}$$

and then obtained (L _ν E ² _ν)_{crit, 2D} ≈ 0.75(L _ν E ² _ν)_{crit, 1D} in rough agreement with simulations using a model for the saturation of non-radial fluid motions (see Section 3.3.3).

Nonetheless, there is likely no monocausal explanation for better heating conditions in multi-D. Yamasaki & Yamada (2006) found, for example, that convective energy transport from the gain radius to the shock also reduces the critical luminosity (although they somewhat overestimated the effect by assuming constant entropy in the entire gain region). Convective energy transport reduces the slope of the pressure gradient between the gain radius (where the pressure is set by the neutrino luminosity and mean energy) and the shock, and thus pushes the shock out by increasing the thermal post-shock pressure. That this effect also plays a role alongside the turbulent pressure can be substantiated by an analysis of neutrino hydrodynamics simulations (Bollig et al. in preparation).

Only a detailed analysis of the properties of turbulence in the gain region (Murphy & Meakin 2011) combined with a model for the interaction of turbulence with a non-spherical accretion shock will reveal the precise combination of multi-D effects that conspire to increase the shock radius compared to 1D. This is no prerequisite for understanding the impact of multi-D effects on the runaway condition as encapsulated by a phenomenological correction factor in Equation (13), since effects like turbulent energy transport, turbulent bulk viscosity, etc. will also scale with the square of the turbulent Mach number in the post-shock region just like the turbulent pressure. They are effectively lumped together in the correction factor (1 + 4/3⟨Ma²⟩)^−3/5. The turbulent Mach number in the post-shock region is thus the crucial parameter for the reduction of the critical luminosity in multi-D, although the coefficient of ⟨Ma²⟩ still needs to be calibrated against multi-D simulations (and maybe different in 2D and 3D).

This does not imply, however, that the energetic requirements for runaway shock expansion in multi-D are fundamentally different from 1D: Runaway still occurs roughly once some material in the gain region first acquires positive total (internal+kinetic+potential) energy e _tot; and the required energy input for this ultimately stems from neutrino heating. ⁸

3.3.3 Saturation of instabilities

What complicates the role of multi-D effects in the neutrino-driven mechanism is that the turbulent Mach number in the gain region itself depends on the heating conditions, which modify the growth rates and saturation properties of convection and the SASI. Considerable progress has been made in recent years in understanding this feedback mechanism and the saturation properties of these two instabilities.

The linear phases of convection and the SASI are now rather well understood. The growth rates for buoyancy-driven convective instability are expected to be of order of the Brunt–Väisälä frequency ω_BV, which can be expressed in terms of P, ρ, c _s, and the local gravitational acceleration g as ⁹

(14) $$\begin{equation} \omega _{\rm BV}^2=g \left(\frac{1}{\rho } \frac{\partial \rho }{\partial r}- \frac{1}{\rho c_{\rm s}^2}\frac{\partial P}{\partial r}\right), \end{equation}$$

which becomes positive in the gain region due to neutrino heating. A first-order estimate yields

(15) $$\begin{equation} \omega _{\rm BV}^2 \sim \frac{G M \dot{Q}_\nu }{4\dot{M} r_{\rm g} ^2 c_{\rm s}^2 \left(r_{\rm sh}-r_{\rm g}\right)} \sim \frac{3\dot{Q}_\nu }{4\dot{M} r_{\rm g} \left(r_{\rm sh}-r_{\rm g}\right)}, \end{equation}$$

using c ² _s ≈ GM/(3r _g) at the gain radius (cp. Müller & Janka 2015). An important subtlety is that advection can stabilise the flow so that ω² _BV > 0 is no longer sufficient for instability unless large seed perturbations in density are already present. Instability instead depends on the more restrictive criterion for the parameter χ (Foglizzo et al. 2006):

(16) $$\begin{equation} \chi = \int \limits _{r_{\rm g}}^{r_{\rm sh}}\frac{\omega _{\rm BV}}{|v_r|}{\rm d}r, \end{equation}$$

with χ ≳ 3 indicating convective instability.

The scaling of the linear growth rate ω_SASI of SASI modes is more complicated, since it involves both the duration τ_cyc of the underlying advective-acoustic cycle as well as a quality factor $\mathcal {Q}$ for the conversion of vorticity and entropy perturbations into acoustic perturbation in the deceleration region below the gain region and the reverse process at the shock (Foglizzo et al. 2006, 2007):

(17) $$\begin{equation} \omega _{\rm SASI} \sim \frac{ \ln |\mathcal {Q}|}{\tau _{\rm cyc}}. \end{equation}$$

For realistic models with strong SASI, one finds $\ln |\mathcal {Q}| \sim 2$ (Scheck et al. 2008; Müller et al. 2012b). SASI growth appears to be suppressed for χ ≳ 3 probably because convection destroys the coherence of the waves involved in the advective-acoustic cycle (Guilet et al. 2010). Interestingly, the demarcation line χ = 3 between the SASI- and convection-dominated regimes is also valid in the non-linear regime if χ is computed from the angle- and time-averaged mean flow (Fernández 2012); and both the SASI and convection appear to drive χ close to this critical value (Fernández 2012).

Both in the SASI-dominated regime and the convection-dominated regime, large growth rates are observed in simulations. It only takes a few tens of milliseconds until the instabilities reach their saturation amplitudes. For this reason, the turbulent Mach number and the beneficial effect of multi-D effects on the heating conditions are typically more sensitive to the saturation mechanism than to initial conditions, so that the onset of shock revival is only subject to modest stochastic variations (Summa et al. 2016). Exceptions apply when the heating conditions vary rapidly, e.g., due to the infall of a shell interface or extreme variations in shock radius (as in the light-bulb models of Cardall & Budiardja 2015), and the runaway condition is only narrowly met or missed (Melson et al. 2015a; Roberts et al. 2016).

The saturation properties of convection were clarified by Murphy et al. (2013), who determined that the volume-integrated neutrino heating rate $\dot{Q}_\nu$ and the convective luminosity L _conv in the gain region roughly balance each other. This can be understood as the result of a self-adjustment process of the accretion flow, whereby a marginally stable, quasi-stationary stratification with χ ≈ 3 is established (Fernández 2012). Müller & Janka (2015) showed that this can be translated into a scaling law that relates the average mass-specific neutrino heating rate $\dot{q}_\nu$ in the gain region to the root mean square average δv of non-radial velocity fluctuations:

(18) $$\begin{equation} \delta v \sim \left[ \dot{q}_\nu (r_{\rm sh}-r_{\rm g}) \right]^{1/3}. \end{equation}$$

That a similar scaling should apply in the SASI-dominated regime is not immediately intuitive. Müller & Janka (2015) in fact tested Equation (18) using a SASI-dominated 2D model and argued that self-adjustment of the flow to χ ≈ 3 will result in the same scaling law as for convection-dominated models. However, models suggest that a different mechanism may be at play in the SASI-dominated regime. Simulations are at least equally compatible with the mechanism proposed by Guilet et al. (2010), who suggested that saturation of the SASI is mediated by parasitic instabilities and occurs once the growth rate of the parasite equals the growth rate of the SASI: Assuming that the Kelvin–Helmholtz instability is the dominant parasite, a simple order-of-magnitude estimate for saturation can be obtained by equating ω_SASI and the average shear rate:

(19) $$\begin{equation} \omega _{\rm SASI} \sim \frac{\delta v}{\Lambda }, \end{equation}$$

where Λ is the effective width of the shear layer. Kazeroni, Guilet, & Foglizzo (2016) find that the Kelvin–Helmholtz instability operates primarily in directions where the shock radius is larger, which suggests Λ = r _{sh, max} − r _g. This results in a scaling law that relates the velocity fluctuations to the average radial velocity $\langle v_\text{r} \rangle$ in the gain region:

(20) $$\begin{equation} \delta v \sim \omega _{\rm SASI} \Lambda \sim \frac{\ln |\mathcal {Q}| (r_{\rm sh,max}-r_{\rm g})}{\tau _{\rm adv}} \sim \ln |\mathcal {Q}| \, |\langle v_\text{r} \rangle |, \end{equation}$$

where we assumed τ_cyc ≈ τ_adv. The quality factor $\mathcal {Q}$ can in principle change significantly with time and between different models. Nonetheless, together with the assumption of a roughly constant quality factor, Equation (20) appears to capture the dynamics of the SASI in 3D quite well for a simulation of an $18\,\text{M}_\odot$ progenitor with the coconut-fmt code (Müller & Janka 2015) as illustrated in Figure 4.

Figure 4. Comparison of the root-mean-square average δv of non-radial velocity component in the gain region (black) with two phenomenological models for the saturation of non-radial instabilities in a SASI-dominated 3D model of an $18\,\text{M}_\odot$ star using the coconut-fmt code. The red curve shows an estimate based on Equation (18), which rests on the assumption of a balance between buoyant driving and turbulent dissipation (Murphy et al. 2013; Müller & Janka 2015). The blue curve shows the prediction of Equation (20), which assumes that saturation is regulated by a balance between the growth rate of the SASI and parasitic Kelvin–Helmholtz instabilities (Guilet, Sato, & Foglizzo 2010). Even though Equation (20) assumes a constant quality factor $|\mathcal {Q}|$ to estimate the SASI growth rate, it appears to provide a good estimate for the dynamics of the model. Interestingly, the saturation models for the SASI- and convection dominated regimes give similar results during later phases even though the mechanism behind the driving instability is completely different.

Equation (18) for the convection-dominated regime and Equation (20) apparently predict turbulent Mach numbers in the same ballpark. This can be understood by expressing $\dot{q}_\nu$ in terms of the accretion efficiency $\eta _{\rm acc} =L_\nu /(GM \dot{M}/r_{\rm g})$ and the heating efficiency $\eta _{\rm heat} =\dot{Q}_\nu /L_\nu$ :

(21) $$\begin{eqnarray} \dot{q}_\nu &=&\frac{\dot{Q}_\nu }{M_{\rm g}}= \eta _{\rm heat} \eta _{\rm acc} \frac{GM \dot{M}}{r_{\rm g} M_{\rm g}}= \eta _{\rm heat} \eta _{\rm acc} \frac{GM}{r_{\rm g} \tau _{\rm adv}} \\ \nonumber &=& \eta _{\rm heat} \eta _{\rm acc} \frac{GM}{r_{\rm sh} \tau _{\rm adv}}\frac{r_{\rm sh}}{r_{\rm gain}}. \end{eqnarray}$$

If we neglect the ratio r _sh/r _g and approximate the average post-shock velocity as $| \langle v_{\rm r} \rangle | \approx \beta ^{-1} \sqrt{GM/r_{\rm sh}}$ (where β is the compression ratio in the shock), we obtain

(22) $$\begin{equation} \dot{q}_\nu \sim \eta _{\rm heat} \eta _{\rm acc} \frac{\beta ^2 |\langle v_\text{r} \rangle |^2 }{\tau _{\rm adv}}, \end{equation}$$

and hence

(23) $$\begin{equation} \delta v \sim (\eta _{\rm heat} \eta _{\rm acc} \beta ^2)^{1/3} |\langle v_\text{r} \rangle |. \end{equation}$$

For plausible values (e.g., η_heat = 0.05η_acc = 2, β = 10), one finds $\delta v \sim 2 |\langle v_\text{r} \rangle |$ , i.e., the turbulent Mach number at saturation is of the same order of magnitude in the convection- and SASI-dominated regimes (where at least $\ln |\mathcal {Q}| \mathord {\sim } 2$ can be reached).

Equations (18) and (20) remain order-of-magnitude estimates; and either of the instabilities may be more efficient at pumping energy into non-radial turbulent motions in the gain region, as suggested by the light-bulb models of Fernández (2015) and Cardall & Budiardja (2015). These authors find that the SASI can lower the critical luminosity in 3D considerably further than convection. Fernández (2010) attributes this to the emergence of the spiral mode of the SASI (Blondin & Mezzacappa 2007; Fernández 2010) in 3D, which can store more non-radial kinetic energy than the SASI sloshing mode in 2D, but this has yet to be borne out by self-consistent neutrino hydrodynamics simulations (see Section 3.4 for further discussion).

3.3.4 Why do models explode more easily in 2D than in 3D?

How can one explain the different behaviour of 2D and 3D models in the light of our current understanding of the interplay between neutrino heating, convection, and the SASI? It seems fair to say that we can presently only offer a heuristic interpretation for the more pessimistic evolution of 3D models.

The most glaring difference between 2D and 3D models (especially in the convection-dominated regime) prior to shock revival lies in the typical scale of the turbulent structures, which are smaller in 3D (Hanke et al. 2012; Couch 2013b; Couch & Ott 2015), whereas the inverse turbulent cascade in 2D (Kraichnan 1967) artificially channels turbulent kinetic energy to large scales. This implies that the effective dissipation length (or also the effective mixing length for energy transport) are smaller in 3D, so that smaller dimensionless coefficients C appear in relations like Equation (18),

(24) $$\begin{equation} \delta v =C \left[ \dot{q}_\nu (r_{\rm sh}-r_{\rm g}) \right]^{1/3}, \end{equation}$$

and the turbulent Mach number will be smaller for a given neutrino heating rate. Indeed, for the $18\,\text{M}_\odot$ model shown in Figure 4, we find

(25) $$\begin{equation} \delta v =0.7 \left[ \dot{q}_\nu (r_{\rm sh}-r_{\rm g}) \right]^{1/3} \end{equation}$$

in 3D rather than what Müller & Janka (2015) inferred from 2D models (admittedly using a different progenitor),

(26) $$\begin{equation} \delta v =\left[ \dot{q}_\nu (r_{\rm sh}-r_{\rm g}) \right]^{1/3}. \end{equation}$$

Following the arguments of Müller & Janka (2015) to infer the correction factor $\left(1+\frac{4\langle {\rm Ma}^2 \rangle }{3}\right)^{-3/5}$ for multi-D effects in Equation (13), one would then expect a considerably larger critical luminosity in 3D, i.e., (L _ν E ² _ν)_{crit, 3D} ≈ 0.85(L _ν E ² _ν)_{crit, 1D} instead of (L _ν E ² _ν)_{crit, 2D} ≈ 0.75(L _ν E ² _ν)_{crit, 1D} in 2D.

Such a large difference in the critical luminosity does not tally with the findings of light-bulb models that show that the critical luminosities in 2D and 3D are still very close to each other. This already indicates that more subtle effects may be at play in 3D that almost compensate the stronger effective dissipation of turbulent motions. The fact that simulations typically show transient phases of stronger shock expansion and more optimistic heating conditions in 3D than in 2D (Hanke et al. 2012; Melson et al. 2015b) also points in this direction.

Furthermore, light-bulb models (Handy, Plewa, & Odrzywołek 2014) and multi-group neutrino hydrodynamics simulations (Melson et al. 2015a; Müller 2015) have demonstrated that favourable 3D effects come into play after shock revival. These works showed that 3D effects can lead to a faster, more robust growth of the explosion energy provided that shock revival can be achieved in the first place.

The favourable 3D effects that are responsible for this may already counterbalance the adverse effect of stronger dissipation in the pre-explosion phase to some extent: Energy leakage from the gain region by the excitation of g-modes is suppressed in 3D because the forward turbulent cascade (Melson et al. 2015a) and (at high Mach number) the more efficient growth of the Kelvin–Helmholtz instability (Müller 2015) brake the downflows before they penetrate the convectively stable cooling layer. Moreover, the non-linear growth of the Rayleigh–Taylor instability is faster for three-dimensional plume-like structures than for 2D structures with planar (Yabe, Hoshino, & Tsuchiya 1991; Hecht et al. 1995; Marinak et al. 1995) or toroidal geometry (as in the context of Rayleigh–Taylor mixing in the stellar envelope during the explosion phase; Kane et al. 2000; Hammer, Janka, & Müller 2010), which might explain why 3D models initially respond more strongly to sudden drops in the accretion rate at shell interfaces and exhibit better heating conditions than their 2D counterparts for brief periods. Finally, the difference in the effective dissipation length in 3D and 2D that is reflected by Equations (25) and (26) may not be universal and depend, e.g., on the heating conditions or the χ-parameter; the results of Fernández (2015) in fact demonstrate that under appropriate circumstances more energy can be stored in non-radial motions in 3D than in 2D in the SASI-dominated regime.

3.4 Outlook: Classical ideas for more robust explosions

The existence of several competing—favourable and unfavourable—effects in 3D first-principle models does not change the fundamental fact that they remain more reluctant to explode than their 2D counterparts. This suggests that some important physical ingredient are still lacking in current simulations. Several avenues towards more robust explosion models have recently been explored. Some of the proposed solutions have a longer pedigree and revisit ideas (rapid rotation in supernova cores, enhanced neutrino luminosities) that have been investigated on and off in supernova theory already before the advent of 3D simulations. The more ‘radical’ solution of invoking strong seed perturbations from convective shell burning to boost non-radial instabilities in the post-shock region will be discussed separately in Section 5.

3.4.1 Rotation and beyond

Nakamura et al. (2014) and Janka et al. (2016) pointed out that rapid progenitor rotation can facilitate explosions in 3D. Janka et al. (2016) ascribed this partly to the reduction of the pre-shock infall velocity due to centrifugal forces, which decreases the ram pressure ahead of the shock. Even more importantly, rotational support also decreases the net binding energy |e _tot| per unit mass in the gain region in their models. They derived an analytic correction factor for the critical luminosity in terms of the average-specific angular momentum j in the infalling shells:

(27) $$\begin{equation} (L_\nu E_\nu ^2)_{\rm crit,rot} \approx (L_\nu E_\nu ^2)_{\rm crit} \times \left(1-\frac{j^2}{2 G M r_{\rm sh}}\right)^{3/5}. \end{equation}$$

Assuming rapid rotation with j ≫ 10¹⁶ cm² s⁻¹, one can obtain a significant reduction of the critical luminosity by several 10% as Janka et al. (2016) tested in a simulation with a modified rotation profile. ¹⁰ For very rapid rotation, other explosion mechanisms also become feasible, such as the magnetorotational mechanism (Akiyama et al. 2003; Burrows et al. 2007b; Winteler et al. 2012; Mösta et al. 2014), or explosions driven by the low-T/W spiral instability (Takiwaki, Kotake, & Suwa 2016).

However, current stellar evolution models do not predict the required rapid rotation rates for these scenarios for the generic progenitors of type IIP supernovae. The typical specific angular momentum at a mass coordinate of $m=1.5\,\text{M}_\odot$ is only of the order of j ~ 10¹⁵ cm² s⁻¹ in models (Heger, Woosley, & Spruit 2005) that include angular momentum transport by magnetic fields generated by the Tayler–Spruit dynamo (Spruit 2002), and asteroseismic measurements of core rotation in evolved low-mass stars suggest that the spin-down of the cores may be even more efficient (Cantiello et al. 2014). For such slow rotation, centrifugal forces are negligible; Equation (27) suggests a change of the critical luminosity on the per-mil level. Neither is rotation expected to affect the character of neutrino-driven convection appreciably because the angular velocity Ω in the gain region is too small. The Rossby number is well above unity:

(28) $$\begin{equation} {\rm Ro}\sim \frac{|v_r|}{ (r_{\rm sh}-r_{\rm g}) \Omega } \sim \frac{r_{\rm s}^2}{\tau _{\rm adv} j} \sim 10, \end{equation}$$

assuming typical values of τ_adv ~ 10 ms and r _sh ~ 100 km.

Magnetic field amplification by a small-scale dynamo or the SASI (Endeve, Cardall, & Budiardja Mezzacappa 2010; Endeve et al. 2012) could also help to facilitate shock revival with magnetic fields acting as a subsidiary to neutrino heating but without directly powering the explosion as in the magnetorototational mechanism. The 2D simulations of Obergaulinger et al. (2014) demonstrated that magnetic fields can help organise the flow into large-scale modes and thereby allow earlier explosions, though the required initial field strengths for this are higher ( $\mathord {\sim } 10^{12} \, {\rm G}$ ) than the typical values predicted by stellar evolution models.

3.4.2 Higher neutrino luminosities and mean energies?

Another possible solution for the problem of missing or delayed explosions in 3D lies in increasing the electron flavour luminosity and mean energy. This is intuitive from Equation (13), where a mere change of $\mathord {\sim }5\%$ in both L _ν and E _ν results in a net effect of 16%, which is almost on par with multi-D effects.

The neutrino luminosity is directly sensitive to the neutrino opacities, which necessitates precision modelling in order to capture shock propagation and heating correctly (Lentz et al. 2012a, 2012b; Müller et al. 2012a; see also Section 4), as well as to other physical ingredients of the CCSN problem that influence the contraction of the proto-neutron star, such as general relativity and the nuclear equation of state (Janka 2012; Müller et al. 2012a; Couch 2013a; Suwa et al. 2013; O’Connor & Couch 2015). Often such changes to the neutrino emission come with counterbalancing side effects (Mazurek’s law); e.g., stronger neutron star contraction will result in higher neutrino luminosities and mean energies, but will also result in a more tightly bound gain region, which necessitates stronger heating to achieve shock revival.

That the lingering uncertainties in the microphysics may nonetheless hold the key to more robust explosions has long been recognised in the case of the equation of state. Melson et al. (2015b) pointed out that missing physics in our treatment of neutrino-matter interactions may equally well be an important part of the solution of the problem shock revival. Exploring corrections to neutral-current scattering cross-section due to the ‘strangeness’ of the nucleon, they found that changes in the neutrino cross-section on the level of a few 10% were sufficient to tilt the balance in favour of explosion for a $20\,\text{M}_\odot$ progenitor. Whilst Melson et al. (2015b) deliberately assumed a larger value for the contribution of strange quarks to the axial form factor of the nucleon than currently measured (Airapetian et al. 2007), the deeper significance of their result is that Mazurek’s law can sometimes be circumvented so that modest changes in the neutrino opacities still exert an appreciable effect on supernova dynamics. A re-investigation of the rates currently employed in the best supernova models for the (more uncertain) neutrino interaction processes that depend strongly on in-medium effects (charged-current absorption/emission, neutral current scattering, Bremsstrahlung; Burrows & Sawyer 1998, 1999; Reddy et al. 1999; Hannestad & Raffelt 1998) may thus be worthwhile (see Bartl, Pethick, & Schwenk 2014; Rrapaj et al. 2015; Shen & Reddy 2014 for some recent efforts).

4 ASSESSMENT OF SIMULATION METHODOLOGY

Considering what has been pointed out in Section 3—the crucial role of hydrodynamic instabilities and the delicate sensitivity of shock revival to the neutrino luminosities and mean energies—it is natural to ask: What are the requirements for modelling the interplay of the different ingredients of the neutrino-driven mechanism accurately? This question is even more pertinent considering that the enormous expansion of the field during the recent years has sometimes produced contradictory results, debates about the relative importance of physical effects, and controversies about the appropriateness of certain simulation methodologies.

Ultimately, only the continuous evolution of the simulation codes, the inclusion of similar physics by different groups, and carefully designed cross-comparisons will eventually produce a ‘concordance model’ of the neutrino-driven mechanism and confirm that simulation results are robust against uncertainties. For 1D neutrino hydrodynamics simulations, this has largely been achieved in the wake of the pioneering comparison paper of Liebendörfer et al. (2005), which has served as reference for subsequent method papers and sensitivity studies in 1D (Müller et al. 2010; Lentz et al. 2012a, 2012b; O’Connor 2015; Just et al. 2015; Summa et al. 2016). Similar results of the Garching-QUB collaboration (Summa et al. 2016) and O’Connor & Couch (2015) with multi-group neutrino transport indicate a trend to a similar convergence in 2D, and more detailed comparisons are underway (see, e.g., https://www.authorea.com/users/1943/articles/97450/_show_article for efforts coordinated by E. O’Connor). Along the road to convergence, it appears useful to provide a preliminary review of some issues concerning the accuracy and reliability of supernova simulations.

4.1 Hydrodynamics

Recently, the discussion of the fidelity of the simulations has strongly focussed on the hydrodynamic side of the problem. As detailed in Section 3, multi-D effects play a crucial role in the explosion mechanism, and are regulated by a balance of driving (by neutrino heating through buoyancy, or by an inherent instability of the flow like the SASI) and dissipation.

4.1.1 Turbulence in supernova simulations

This balance needs to be modelled with sufficient physical and numerical accuracy. On the numerical side, the challenge consists in the turbulent high-Reynolds number flow, and the question arises to what extent simulations with relatively coarse resolution can capture this turbulent flow accurately. Various authors (Handy et al. 2014; Abdikamalov et al. 2015; Radice, Couch, & Ott 2015; Roberts et al. 2016) have stressed that the regime of fully developed turbulence cannot be reached with the limited resolution affordable to cover the gain region ( $\mathord {\sim } 100$ zones, or even less) in typical models, and Handy et al. (2014) thus prefer to speak of ‘perturbed laminar flow’ in simulations. Attempts to quantify the effective Reynolds number of the flow using velocity structure functions and spectral properties of the post-shock turbulence (Handy et al. 2014; Abdikamalov et al. 2015; Radice et al. 2015) put it at a few hundred at best, and sometimes even below 100.

This is in line with rule-of-thumb estimates based on the numerical diffusivity for the highest wavenumber (odd-even) modes in Godunov-based schemes as used in many supernova codes. This diffusivity can be calculated analytically (Appendix D of Müller 2009; see also Arnett & Meakin 2016 for a simpler estimate). For Riemann solvers that take all the wave families into account (e.g., Colella & Glaz 1985; Toro, Spruce, & Speares 1994; Mignone & Bodo 2005; Donat & Marquina 1996), the numerical kinematic viscosity ν_num in the subsonic regime is roughly given in terms of the typical velocity jump per cell δv _gs and the cell width δl as ν_num ~ δl δv _gs. Relating δv _gs to the turbulent velocity v and scale l of the largest eddy as δv _gs ~ v(δl/l)^1/3 (i.e., assuming Kolmogorov scaling) yields a numerical Reynolds number of

(29) $$\begin{equation} {\rm Re} = \frac{v l}{\nu _{\rm num}} \sim \left(\frac{l}{\delta l}\right)^{4/3}=N^{4/3}, \end{equation}$$

where N is the number of zones covering the largest eddy scale. For more diffusive solvers like HLLE (Einfeldt 1988), one obtains ν_num ~ δl c _s ~ δl v Ma⁻¹ instead and

(30) $$\begin{equation} {\rm Re} \sim (l/\delta l) {\rm Ma} \sim N \, {\rm Ma}, \end{equation}$$

i.e., such solvers are strongly inferior for subsonic flow with low Mach number Ma.

Such coarse estimates are to be taken with caution, however. The numerical dissipation is non-linear and self-regulated as typical of implicit large-eddy simulations (ILES, Boris et al. 1992; Grinstein, Margolin, & Rider 2007). In fact, the estimates already demonstrate that simply comparing the resolution in codes with different solvers and grid geometries can be misleading. Codes with three-wave solvers like Vertex-Prometheus (Rampp & Janka 2002; Buras et al. 2006a) and coconut-fmt (Müller & Janka 2015) of the MPA-QUB collaboration, flash (Fryxell et al. 2000) as used in Couch (2013a) and subsequent work by S. Couch and E. O’Connor, and the vh-1 hydro module (Blondin, Stevens, & Kallman 1991) in the Chimera code of the Oak Ridge-Florida Atlantic-NC State collaboration, have less stringent resolution requirements than HLLE-based codes (Ott et al. 2012; Kuroda, Kotake, & Takiwaki 2012). The reconstruction method, special tweaks for hydrostatic equilibrium (or an the lack of such a treatment), as well as the grid geometry and grid-induced perturbations (Janka et al. 2016; Roberts et al. 2016) also affect the behaviour and resolution-dependence of the simulated turbulence.

4.1.2 Resolution requirements—a critical assessment

Regardless of the employed numerical schemes, the fact remains that the achievable numerical Reynolds number in supernova simulations is limited, and that the regime of fully developed turbulence (Re ≫ 1000) will not be achieved in the near future, as it would require ≳ 512 radial zones in the gain region alone. The question for supernova models, however, is not whether all the facets of turbulence in inviscid flow can be reproduced, but whether the flow properties that matter for the neutrino-driven mechanism are computed with sufficient accuracy. In fact, one cannot even hope that simply cranking up the numerical resolution with ILES methods would give the correct solution: In reality, non-ideal effects such as neutrino viscosity and drag (van den Horn & van Weert 1984; Burrows 1988; Jedamzik, Katalinić, & Olinto 1998; Guilet, Müller, & Janka 2015) come into play, and deviations of the turbulent Prandtl number from unity as well as MHD effects like a small-scale dynamo (see Section 3.4) can complicate the picture even for non-rotating, weakly magnetised supernova cores. These effects will likely not grossly alter the dynamics of convection and the SASI, but the physical reality may be slightly different from the limit of infinite resolution if these effects are not accounted for and inviscid flow is assumed instead.

At the end of the day, these additional complications and the finite resolution probably have a limited effect on supernova dynamics, since they only affect a correction term to the critical luminosity such as (1 + 4/3⟨Ma²⟩)^−3/5 in Equation (13) through the effective dissipation length that determines the non-dimensional coefficient in Equation (18). If we repeat the analytic estimate for L _crit of Müller & Janka (2015), but assume stronger dissipation and decrease their critical Mach number at shock revival Ma² _crit = 0.4649 by 10%, then Equation (13) suggests an increase of the critical luminosity from 74.9% of the 1D value to of 76.6% of the 1D value, which is a minute change. Modelling turbulent dissipation within 10% uncertainty thus seems wholly sufficient given that one can hardly hope to achieve 1% accuracy for the neutrino luminosities and mean energies.

The turbulent dissipation does not change without bounds with increasing resolution, but eventually reaches an asymptotic limit at high Reynolds numbers. Although most supernova simulation may not fully reach this asymptotic regime, they do not fall far short of it: The works of Handy et al. (2014) and Radice et al. (2015, 2016) suggest that this level of accuracy in the turbulent dissipation can be reached even with moderate resolution (< 100 grid points per direction, ~ 2° resolution in angle in spherical polar coordinates) in the gain region with higher order reconstruction methods and accurate Riemann solvers. Problems due to stringent resolution requirements may still lurk elsewhere, though, e.g., concerning SASI growth rates as already pointed out 10 yr ago by Sato, Foglizzo, & Fromang (2009). Resolution studies and cross-comparisons thus remain useful, though cross-comparisons are of course hampered by the different physical assumptions used in different codes and the feedback processes in the supernova core. For this reason, a direct comparison of, e.g., turbulent kinetic energies and Mach numbers between different models is not necessarily meaningful. The dimensionless coefficients governing the dynamics of non-radial instabilities such the proportionality constant $\eta _{\rm conv}=v_{\rm turb}/[\dot{q_\nu } (r_{\rm sh}-r_{\rm g})]$ in Equation (18) or the quality factor $\mathcal {Q}$ in Equation (17) may be more useful metrics of comparison.

4.2 Neutrino transport

The requirements on the treatment of neutrino heating and cooling are highly problem-dependent. The physical principles behind convection and the SASI can be studied with simple heating and cooling functions in a light-bulb approach, and such an approach is indeed often advantageous as it removes some of the feedback processes that complicate the analysis of full-scale supernova simulations. To model the fate and explosion properties of concrete progenitors in a predictive manner, some form of neutrino transport is required, and depending on the targeted level of accuracy, the requirements become more stringent; e.g., higher standards apply when it comes to predicting supernova nucleosynthesis. There is no perfect method for neutrino transport in supernovae as yet. Efforts towards a solution of the full 6D Boltzmann equation are underway (e.g., Cardall, Endeve, & Mezzacappa 2013; Peres et al. 2014; Radice et al. 2013; Nagakura, Sumiyoshi, & Yamada 2014), but not yet ripe for real supernova simulations.

Neutrino transport algorithms (beyond fully parameterised light-bulb models) currently in use for 1D and multi-D models include:

• • leakage schemes as, e.g., in O’Connor & Ott (2010, 2011), Ott et al. (2013), and Couch & O’Connor (2014);

• • the isotropic diffusion source approximation (IDSA) of Liebendörfer, Whitehouse, & Fischer (2009);

• • one-moment closure schemes employing prescribed flux factors (Scheck et al. 2006), flux-limited diffusion as in the Vulcan code (Livne et al. 2004; Walder et al. 2005), the Chimera code (Bruenn 1985; Bruenn et al. 2013), and the Castro code (Zhang et al. 2013; Dolence, Burrows, & Zhang 2015), or a dynamic closure as in the coconut-fmt code;

• • two-moment methods employing algebraic closures in 1D (O’Connor 2015) and multi-D (Obergaulinger & Janka 2011; Kuroda et al. 2012; Just et al. 2015; Skinner et al. 2016; O’Connor & Couch 2015; Roberts et al. 2016; Kuroda, Takiwaki, & Kotake 2016) or variable Eddington factors from a model Boltzmann equation (Burrows et al. 2000b; Rampp & Janka 2002; Buras et al. 2006a; Müller et al. 2010);

• • discrete ordinate methods for the Boltzmann equation, mostly in 1D (Mezzacappa & Bruenn 1993; Yamada, Janka, & Suzuki 1999; Liebendörfer et al. 2004) or, at the expense of other simplifications, in multi-D (Livne et al. 2004; Ott et al. 2008; Nagakura et al. 2016; only for static configurations: Sumiyoshi et al. 2015).

This list should not be taken as a hierarchy of accuracy; it mere reflects crudely the rigour in treating one aspect of the neutrino transport problem, i.e., the angle-dependence of the radiation field in phase space. When assessing neutrino transport methodologies, there are other, equally important factors that need to be taken into account when comparing different modelling approaches.

Most importantly, the sophistication of the microphysics varies drastically. On the level of one-moment and two-moment closure models, it is rather the neutrino microphysics that decides about the quantitative accuracy. The 3D models of the MPA-QUB group (Melson et al. 2015a, 2015b; Janka et al. 2016) and the Chimera team (Lentz et al. 2015) currently represent the state-of-the-art in this respect; though other codes (O’Connor 2015; Just et al. 2015; Skinner et al. 2016; Kuroda et al. 2016) come close.

Often, the neutrino physics is simplified considerably, however. Some simulations disregard heavy flavour neutrinos altogether (e.g., Suwa et al. 2010; Takiwaki et al. 2012), or only treat them by means of a leakage scheme (Takiwaki et al. 2014; Pan et al. 2016). This affects the contraction of the proto-neutron star and thus indirectly alters the emission of electron flavour neutrinos and the effective inner boundary for the gain region as well.

Amongst multi-D codes, energy transfer due to inelastic neutrino-electron scattering (NES) is routinely taken into account only in the Vertex code (Rampp & Janka 2002; Buras et al. 2006a; Müller et al. 2010) of the MPA-QUB collaboration, the Alcar code (Just et al. 2015), the Chimera code of the Chimera team (Bruenn 1985; Bruenn et al. 2013), and the Fornax code of the Princeton group (Skinner et al. 2016). Without NES (Bruenn 1985) and modern electron capture rates (Langanke et al. 2003), the core mass at bounce is larger and the shock propagates faster at early times (Lentz et al. 2012a, 2012b). In multi-D, this can lead to unduly strong prompt convection. Because of this problem, a closer look at the bounce dynamics is in order whenever explosions occur suspiciously early (< 100 ms after bounce). Parameterising deleptonisation during collapse (Liebendörfer 2005) provides a workaround to some extent.

The recoil energy transfer in neutrino-nucleon scattering effectively reshuffles heavy flavour neutrino luminosity to electron flavour luminosity in the cooling region (Müller et al. 2012a) and hence critically influences the heating conditions in the gain region. Amongst multi-D codes, only Vertex and Chimera currently take this into account, and the code coconut-fmt (Müller & Janka 2015) uses an effective absorption opacity for heavy flavour neutrinos to mimic this phenomenon.

Vertex and Chimera are also the only multi-D codes to include the effect of nucleon–nucleon correlations (Burrows & Sawyer 1998, 1999; Reddy et al. 1999) on absorption and scattering opacities. Nucleon correlations have a huge impact during the cooling phase, which they shorten by a factor of several (Hüdepohl et al. 2009). Their role during the first second after bounce is not well explored. Considering that the explosion energetics are determined on a timescale of seconds (Müller 2015; Bruenn et al. 2016), it is plausible that the increased diffusion luminosity from the neutron star due to in-medium corrections to the opacities may influence the explosion energy to some extent.

Gray schemes (Fryer & Warren 2002; Scheck et al. 2006; Kuroda et al. 2012) cannot model neutrino heating and cooling accurately; an energy-dependent treatment is needed because of the emerging neutrino spectra are highly non-thermal with a pinched high-energy tail (Janka & Hillebrandt 1989; Keil, Raffelt, & Janka 2003).

Some multi-D codes use the ray-by-ray-plus approximation (Buras et al. 2006a), which exaggerates angular variations in the radiation field, and has been claimed to lead to spuriously early explosions in some cases in conjunction with artificially strong sloshing motions in 2D (Skinner et al. 2016). Whether this is a serious problem is unclear in the light of similar results of Summa et al. (2016) for ray-by-ray-plus models and O’Connor & Couch (2015) for fully two-dimensional two-moment transport. On the other hand, fully multi-dimensional flux limited diffusion approaches smear out angular variations in the radiation field too strongly (Ott et al. 2008).

Neglecting all or part of the velocity-dependent terms in the transport equations potentially has serious repercussions. Neglecting only observer correction (Doppler shift, compression work, etc.) as, e.g., in Livne et al. (2004) can already have an appreciable impact on the dynamics (Buras et al. 2006a; Lentz et al. 2012a). Disregarding even the co-advection of neutrinos with the fluid (O’Connor 2015; Roberts et al. 2016) formally violates the diffusion limit and effectively results in an extra source term in the optically thick regime due to the equilibration of matter with lagging neutrinos:

(31) $$\begin{equation} \dot{q}_\nu \approx \rho ^{-1} \mathbf {v} \cdot \nabla E_{\rm eq}, \end{equation}$$

where E _eq is the equilibrium neutrino energy density. Judging from the results of O’Connor & Couch (2015) and Roberts et al. (2016), which are well in line with results obtained with other codes, the effect may not be too serious in practice, though. It should also be noted that (semi-)stationary approximations of the transport equation (Liebendörfer et al. 2009; Müller & Janka 2015) avoid this problem even if advection terms are not explicitly included.

Leakage-based schemes as used, e.g., in Ott et al. (2012), Couch & Ott (2015), Abdikamalov et al. (2015), and Couch et al. (2015) also manifestly fail to reproduce the diffusion limit. Here, however, the violation of the diffusion limit is unmistakable and can severely affect the stratification of the gain region and, in particular, the cooling region. Together with ad hoc choices for the flux factor for calculating the heating rate, this can result in inordinately high heating efficiencies immediately after bounce and a completely inverted hierarchy of neutrino mean energies. It compromises the dynamics of leakage models to an extent that they can only be used for very qualitative studies of the multi-D flow in the supernova core.

There is in fact no easy lesson to be learned from the pitfalls and complications that we have outlined. In many contexts, approximations for the neutrino transport are perfectly justified for a well-circumscribed problem, and feedback processes sometimes mitigate the effects of simplifying assumptions. It is crucial, though, to be aware of the impact that such approximations can potentially have, and our (incomplete) enumeration is meant to provide some guidance in this respect.

5 FUTURE DIRECTIONS: MULTI-D EFFECTS IN SUPERNOVA PROGENITORS

Given the sophisticated simulation methodology employed in the best currently available supernova codes, one may be tempted to ask whether another missing ingredient for robust neutrino-driven explosion is to be sought elsewhere. One recent idea, first proposed by Couch & Ott (2013), focusses on the progenitor models used in supernova simulations. The twist consists in an extra ‘forcing’ of the non-radial motions in the gain region by large seed perturbations in the infalling shells. Such seed perturbations will arise naturally in active convective burning shells (O burning, and perhaps also Si burning) that reach the shock during the first few hundred milliseconds after bounce.

5.1 Role of pre-collapse perturbations in the neutrino-driven mechanism

In default of multi-D progenitor models, this new variation of the neutrino-driven mechanism was initially studied by imposing large initial perturbations by hand in leakage-based simulations (Couch & Ott 2013, 2015) and multi-group neutrino hydrodynamics simulations (Müller & Janka 2015); the earlier light-bulb-based models of Fernández (2012) also touched parts of the problem. The results of these investigations were mixed, even though some of these calculations employed perturbations far in excess of what estimates based on mixing-length theory (Biermann 1932; Böhm-Vitense 1958) suggest: For example, Couch & Ott (2013) used transverse velocity perturbations with a peak Mach number of Ma = 0.2 in their 3D models, and found a small beneficial effect on shock revival, which, however, was tantamount to a change of the critical neutrino luminosity by only $\mathord {\sim } 2\%$ . The more extensive 2D parameter study of different solenoidal and compressive velocity perturbations and density perturbations by Müller & Janka (2015) established that both significant perturbation velocities (Ma ≳ 0.1) as well as large-scale angular structures (angular wavenumber ℓ ≲ 4) need to be present in active convective shell in order to reduce the critical luminosity appreciably, i.e., by ≳ 10%.

These parametric studies already elucidated the physical mechanism whereby pre-collapse perturbations can facilitate shock revival. Müller & Janka (2015) highlighted the importance both of the infall phase as well as the interaction of the perturbations with the shock. Linear perturbation theory shows that the initial perturbations are amplified during collapse (Lai & Goldreich 2000; Takahashi & Yamada 2014). This not only involves a strong growth of transverse velocity perturbations as δv_t ∝r ⁻¹, but even more importantly a conversion of the initially dominating solenoidal velocity perturbations with Mach number Ma_conv into density perturbations δρ/ρ ≈ Ma (Müller & Janka 2015) during collapse, i.e., the relative density perturbations are much larger ahead of the shock than during quasi-stationary convection, where δρ/ρ ≈ Ma². ¹¹

Large density perturbations ahead of the shock imply a pronounced asymmetry in the pre-shock ram pressure and deform the shock, creating fast lateral flows as well as post-shock density and entropy perturbations that buoyancy then converts into turbulent kinetic energy. The direct injection of kinetic energy due to infalling turbulent motions may also play a role (Abdikamalov et al. 2016), though it appears to be subdominant (Müller & Janka 2015; Müller et al. 2016a). A very crude estimate for the generation of additional turbulent kinetic energy due to the different processes as well as turbulent damping in the post-shock region has been used by Müller et al. (2016a) to estimate the reduction of the critical luminosity as

(32) $$\begin{equation} (L_\nu E_\nu ^2)_{\rm crit,pert} \approx (L_\nu E_\nu ^2)_{\rm crit,3D} \left( 1- 0.47 \frac{{\rm Ma_{\rm conv}}}{\ell \eta _{\rm acc} \eta _{\rm heat}} \right), \end{equation}$$

in terms of the pre-collapse Mach number Ma_conv of eddies from shell burning, their typical angular wavenumber ℓ, and the accretion efficiency $\eta _{\rm acc}=L_\nu /(GM \dot{M} r_{\rm gain})$ and heating efficiency η_heat during the pre-explosion phase.

A more rigorous understanding of the interaction between infalling perturbations, the shock, and non-radial motions in the post-shock region is currently emerging: Abdikamalov et al. (2016) studied the effect of upstream perturbations on the shock using the linear interaction approximation of Ribner (1953) and argue, in line with Müller et al. (2016a), that a reduction of the critical luminosity by > 10% is plausible. Their estimate may, however, be even too pessimistic as they neglect acoustic perturbations upstream of the shock. Different from Abdikamalov et al. (2016), the recent analysis of Takahashi et al. (2016) also takes into account that instabilities or stabilisation mechanisms operate in the post-shock flow, and studied the (linear) response of convective and SASI eigenmodes to forcing by infalling perturbations. A rigorous treatment along these lines that explains the saturation of convective and SASI modes as forced oscillators with non-linear damping remains desirable.

5.2 The advent of 3D supernova progenitor models

The parametric studies of Couch & Ott (2013, 2015) and Müller & Janka (2015) still hinged on uncertain assumptions about the magnitude and scale of the seed perturbations left by O and Si shell burning. Various pioneering studies of advanced shell burning stages (O, Si, C burning) (Arnett 1994; Bazan & Arnett 1994, 1998; Asida & Arnett 2000; Kuhlen, Woosley, & Glatzmaier 2003; Meakin & Arnett 2006, 2007b, 2007a; Arnett & Meakin 2011; Viallet et al. 2013; Chatzopoulos, Graziani, & Couch 2014) merely indicated that convective Mach numbers of a few 10⁻² and the formation of large-scale eddies are plausible, but did not permit a clear-cut judgement about whether pre-collapse perturbations play a dynamical role in the neutrino-driven mechanism.

The situation has changed recently with the advent of models of convective shell burning that have been evolved up to collapse. The idea here is to calculate the last few minutes prior to collapse to obtain multi-dimensional initial conditions, whilst ignoring potential long-term effects in 3D such as convective boundary mixing (which we discuss in Section 5.3). Couch et al. (2015) performed a 3D simulation of the last minutes of Si shell burning in a $15\,\text{M}_\odot$ star. The simulation was limited to an octant, and nuclear quasi-equilibrium during Si burning was only treated with a small network. More importantly, the evolution towards collapse was artificially accelerated by artificially increasing electron capture rates in the iron core. As pointed out by Müller et al. (2016a), this can alter the shell evolution and the convective velocities considerably. Since the shell configuration and structure at collapse varies considerably in 1D models, such an exploratory approach is nonetheless still justified (see below).

Müller et al. (2016a) explored the more generic case where Si shell burning is extinguished before collapse and the O shell is the innermost active convective region. In their 3D simulation of the last 5 min of O shell burning in an $18\,\text{M}_\odot$ progenitor, they circumvented the aforementioned problems by excising the non-convective Fe and Si core and contracting it in accordance with a 1D stellar evolution model. Moreover, Müller et al. (2016a) simulated the entire sphere using an overset Yin–Yang grid (Kageyama & Sato 2004; Wongwathanarat, Hammer, & Müller 2010) as implemented (with some improvements) in the Prometheus supernova code (Melson 2013; Melson et al. 2015a).

The implications of these simulations for supernova modelling are mixed. The typical convective Mach number in Couch et al. (2015) was only $\mathord {\sim }0.02$ , and whilst they found large-scale motions, the scale of the pre-collapse perturbations was still limited by the restriction to octant symmetry. Perturbations of such a magnitude are unlikely to reduce the critical luminosity considerably (Section 5.1). Consequently, supernova simulations starting from 1D and 3D initial conditions using a leakage scheme performed by Couch et al. (2015) did not show a qualitative difference; both 1D and 3D initial conditions result in explosions, though the shock expands slightly faster in the latter case. The use of a leakage scheme and possible effects of stochasticity preclude definite conclusions from these first results.

The typical convective Mach number in the $18\,\text{M}_\odot$ model of Müller et al. (2016a) is considerably larger ( $\mathord {\sim }0.1$ ), and their simulation also showed the emergence of a bipolar (ℓ = 2) flow structure, which lead them to predict a relatively large reduction of the critical luminosity by 12. . .24%, which would accord a decisive role to 3D initial conditions in the neutrino-driven mechanism at least in some progenitors. A first 3D multi-group neutrino hydrodynamics simulation of their $18\,\text{M}_\odot$ progenitor using the coconut-fmt code appears to bear this out (Müller et al. in preparation): Figure 5 shows the shock radius both for two simulations using 3D and 1D initial conditions, respectively: In the former case, shock revival occurs around 250 ms after bounce thanks to the infall of the convectively perturbed oxygen shell, whereas no explosion develops in the reference simulation by the end of the run more than 600 ms after bounce. An analysis of the heating conditions indicates that the non-exploding reference model is clearly not a near miss at 250 ms. The effect of 3D initial conditions is thus unambiguously large and sufficient to change the evolution qualitatively. Moreover, the model indicates that realistic supernova explosion energies are within reach in 3D as well: The diagnostic explosion energy reaches 5 × 10⁵⁰ erg and still continues to mount by the end of the simulation 1.43 s after bounce. It is also interesting to note that the initial asymmetries are clearly reflected in the explosion geometry (Figure 6) as speculated by Arnett & Meakin (2011). Incidentally, the model also shows that the accretion of convective regions does not lead to the formation of the ‘accretion belts’ proposed by Gilkis & Soker (2014) as an ingredient for their jittering-jet mechanism.

Figure 5. Impact of pre-collapse asphericities on shock revival in 3D multi-group neutrino hydrodynamics simulations of an $18\,\text{M}_\odot$ progenitor. The plot shows the minimum, maximum (solid lines), and average (dashed) shock radii for a model using 3D initial conditions (black) from the O shell burning simulation of Müller et al. (2016a) and a spherically averaged version of the same progenitor (red). The gain radius (dash-dotted) and the proto-neutron star radius (dotted, defined by a fiducial density of 10¹¹ g cm⁻³) are shown only for the model starting from 3D initial conditions; they are virtually identical for both models. A neutrino-driven explosion is triggered roughly 0.25 s after bounce aided by the infall of the convectively perturbed oxygen shell in the model using 3D initial conditions. The simulation starting from the 1D progenitor model exhibits steady and strong SASI oscillations after 0.25 s, but does not explode at least for another 0.3 s.

Figure 6. Top row: Radial velocity in units of cm s⁻¹ (top left) and mass fraction of Si (top right) at the onset of collapse in the 3D progenitor model of an $18\,\text{M}_\odot$ star of Müller et al. (2016a). Bottom row: Entropy in units of k_b /nucleon (bottom left) and mass fraction of Si (bottom right) in the ensuing neutrino-driven explosion 1.43 s after bounce from. All plots show equatorial slices from the 3D simulation. It can be seen that the geometry of the initial conditions is still imprinted on the explosion to some extent with stronger shock expansion in the direction of updrafts of Si rich ashes in the O burning shell. This is a consequence of the forced deformation of the shock around the onset of the explosion.

Whether 3D initial conditions generally play an important role in the neutrino-driven mechanism cannot be answered by studying just two progenitors, aside from the fact that the models of Couch et al. (2015) and Müller et al. (2016a) still suffer from limitations. The properties (width, nuclear energy generation rate) and the configuration of convective burning shells at collapse varies tremendously across different progenitors in 1D stellar evolution models as, e.g., the Kippenhahn diagrams in the literature indicate (Heger, Langer, & Woosley 2000; Chieffi & Limongi 2013; Sukhbold & Woosley 2014; Cristini et al. 2016). The interplay of convective burning, neutrino cooling, and the contraction/re-expansion of the core and the shells sometimes leave inversions in the temperature stratification and a complicating layering of material at different nuclear processing stages. For this reason, 1D stellar evolution models sometimes show a highly dynamic behaviour immediately prior to collapse with shells of incompletely burnt material flaring up below the innermost active shell. This is illustrated by follow-up work to Müller et al. (2016a) shown in Figure 7, where a partially processed layer with unburnt O becomes convective shortly before collapse due to violent burning and is about to merge with the overlying O/Ne shell before collapse intervenes.

Figure 7. Radial velocity in units of cm s⁻¹ (shown in 90°-wedges in the left half of each plot) and mass fraction X _O of oxygen during the last minutes of shell burning in an $12.5\,\text{M}_\odot$ progenitor. Snapshots at 175 s (top left), 66 s (top right), 24 s (bottom left) before collapse, and at the onset of collapse (bottom right) are shown. The residual oxygen in a thin, almost O-depleted shell (red) starts to burn vigorously due to the contraction of the core (top right). As the entropy of this shell increases and matches that of an almost unprocessed, O-rich shell (blue) and the active Ne shell (cyan), it expands outwards by ‘encroachment’ (bottom left), but there is insufficient time for the shells to merge completely before collapse (bottom right). Note that this is not a qualitatively new phenomenon in 3D; similar events occur in 1D stellar evolution models.

The diverse shell configurations in supernova progenitors need to be thoroughly explored in 3D before a general verdict on the efficacy of convective seed perturbations in aiding shock revival can be given. Since the bulk properties of the flow (typical velocity, eddy scales) in the interior of the convective shells are apparently well captured by mixing-length theory (Arnett, Meakin, & Young 2009; Müller et al. 2016a), the convective Mach numbers and eddy scales predicted from 1D stellar evolution models can provide guidance for exploring interesting spots in parameter space.

5.3 Convective boundary mixing—how uncertain is the structure of supernova progenitors?

In what we discussed so far, we have considered multi-D effects in advanced convective burning stages merely because of their role in determining the initial conditions for stellar collapse. They could also have an important effect on the secular evolution of massive stars long before the supernova explosion, and thereby change critical structural properties of the progenitors, such as the compactness parameter (O’Connor & Ott 2011). Whilst mixing-length theory (Biermann 1932; Böhm-Vitense 1958) may adequately describe the mixing in the interior of convective zones, ¹² the mixing across convective boundaries is less well understood, and may play an important role in determining the pre-collapse structure of massive stars along with other non-convective processes (e.g., Heger et al. 2000; Maeder & Meynet 2004; Heger et al. 2005; Young et al. 2005; Talon & Charbonnel 2005; Cantiello et al. 2014) for mixing and angular momentum transport. That some mixing beyond the formally unstable regions needs to be included has long been known (Kippenhahn & Weigert 1990). Phenomenological recipes for this include extending the mixed region by a fraction of the local pressure scale height, or adding diffusive mixing in the formally stable regions with a calibrated functional dependence on the distance to the boundary (Freytag, Ludwig, & Steffen 1996; Herwig et al. 1997).

The dominant mechanism for convective boundary mixing during advanced burning stages is entrainment (Fernando 1991; Meakin & Arnett 2007b; Viallet et al. 2015) due to the growth of the Kelvin–Helmholtz or Holmböe instability at the shell interfaces. For interfaces with a discontinuous density jump as often encountered in the interiors of evolved massive stars, the relevant dimensionless number for such shear-driven instabilities is the bulk Richardson number Ri_B. For entrainment driven by turbulent convection, one has

(33) $$\begin{equation} {\rm Ri}_{\rm B}=\frac{g l\, \delta \rho /\rho }{v_{\rm conv}^2}, \end{equation}$$

in terms of the local gravitational acceleration g, the density contrast δρ/ρ at the interface, the typical convective velocity v _conv in the convective region, and the integral scale l of the convective eddies. Equating l with the pressure scale height l = P/ρg allows us to re-express Ri_B in terms of the convective Mach number Ma_conv and the adiabatic exponent γ:

(34) $$\begin{equation} {\rm Ri}_{\rm B}= \frac{\delta \rho }{\rho } \frac{g l}{v_{\rm conv}^2} = \frac{\delta \rho }{\rho } \frac{P}{\rho v_{\rm conv}^2} = \frac{\delta \rho }{\rho } \frac{1}{\gamma {\rm Ma}_{\rm conv}^2}. \end{equation}$$

Deep in the stellar core, Ma_conv is typically small during most evolutionary phases, and Ri_B is large so that the convective boundaries are usually very ‘stiff’ (Cristini et al. 2016).

Various power laws for the entrainment rate have been proposed in the general fluid dynamics literature (Fernando 1991; Strang & Fernando 2001) and astrophysical studies (Meakin & Arnett 2007b) of interfacial mixing driven by turbulent convection on one side of the interface. In the astrophysical context, it is convenient to translate these into a power law for the mass flux $\dot{M}_{\rm entr}$ of entrained material into the convective region,

(35) $$\begin{equation} \dot{M}_{\rm entr} =4 \pi r^2 \rho v_{\rm conv} A \, {\rm Ri}_{\rm B}^{-n}, \end{equation}$$

with a proportionality constant A and a power-law exponent n. Here, ρ is the density on the convective side of the interface.

A number of laboratory studies (Fernando 1991; Strang & Fernando 2001) and astrophysical simulations (Meakin & Arnett 2007b; Müller et al. 2016a) suggest values of A ~ 0.1 and n = 1. This can be understood heuristically by assuming that layer of width δl ~ Av ² _conv/(g δρ/ρ) always remains well mixed, ¹³ and that a fraction δl/l of the mass flux $\dot{M}_{\rm down} =2 \pi r^2 \rho v_{\rm conv}$ in the convective downdrafts comes from this mixed layer.

This estimate is essentially equivalent to another one proposed in a slightly different context (ingestion of unburnt He during core-He burning; Constantino et al. 2015) by Spruit (2015), who related the ingestion (or entrainment) rate into a convective zone to the convective luminosity L _conv. Spruit’s argument can be interpreted as one based on energy conservation; work is needed to pull material with positive buoyancy from an outer shell down into a deeper one, and the energy that is tapped for this purpose comes from convective motions. Since L _conv ~ 4πr ²ρv ³ _conv, we can write Equation (35) as

(36) $$\begin{equation} \dot{M}_{\rm entr} =A \times \frac{4 \pi r^2 \rho v_{\rm conv}^3 }{g l\,\delta \rho /\rho } \approx A \times \frac{L_{\rm conv}}{g l \, \delta \rho /\rho }, \end{equation}$$

which directly relates the entrainment rate to the ratio of L _conv and the potential energy of material with positive buoyancy after downward mixing over an eddy scale l. The entrainment law (35), the argument of Spruit (2015), and the proportionality of the entrainment rate with L _conv found in the recent work of Jones et al. (2016b) on entrainment in highly resolved idealised 3D simulation of O shell burning appear to be different sides of the same coin.

5.4 Long-term effects of entrainment on the shell structure?

How much will entrainment affect the shell structure of massive stars in the long term? First numerical experiments based on the entrainment law of Meakin & Arnett (2007b) were performed by Staritsin (2013) for massive stars on the main sequence ¹⁴ and did not reveal dramatic differences in the size of the convective cores compared to more familiar, calibrated recipes for core overshooting.

Taking Equation (36) at face value allows some interesting speculations about the situation during advanced burning stages. Since the convective motions ultimately feed on the energy generated by nuclear burning E _burn, we can formulate a time-integrated version of Equation (36) for the entrained mass ΔM _entr over the life time of a convective shell:

(37) $$\begin{eqnarray} \frac{GM}{r} \frac{\delta \rho }{\rho } \Delta M_{\rm entr} &\lesssim & A E_{\rm burn}, \end{eqnarray}$$

(38) $$\begin{eqnarray} \frac{GM}{r} \frac{\delta \rho }{\rho } \Delta M_{\rm entr} &\lesssim & A M_{\rm shell} \Delta Q, \end{eqnarray}$$

where M _shell is the (final) mass of the shell, and ΔQ is the nuclear energy release per unit mass. With GM/r ~ 2e _int in stellar interiors, we can estimate ΔM _entr in terms ΔQ and the internal energy e _int at which the burning occurs, ¹⁵

(39) $$\begin{equation} \Delta M_{\rm entr} \lesssim A M_{\rm shell} \left(\frac{\delta \rho }{\rho }\right)^{-1} \frac{\Delta Q}{2 e_{\rm int} }. \end{equation}$$

For O burning at ~ 2 × 10⁹ K and with ΔQ ≈ 0.5 MeV/nucleon, the factor ΔQ/(2e _int) is of order unity. Typically, the density contrast δρ/ρ between adjacent shells is also not too far below unity. Since A ≈ 0.1, this suggests that the shell growth due to entrainment comes up to at most a few tens of percent during O shell burning unless δρ/ρ is rather small to begin with. Thus, a result of entrainment might be that convective zones may swallow thin, unburnt shells with a small density contrast before bounce, whereas the large entropy jumps between the major shells are maintained and even enhanced as a result of this cannibalisation.

For C burning, the long-term effect of entrainment could be somewhat larger than for O burning due to the lower temperature threshold and the higher ratio ΔQ/2e _int; for Si burning, the effect should be smaller. During earlier phases, our estimates break down because the convective flux carries only a small fraction of the energy generation by nuclear burning. If this is taken into account, the additional growth of convective regions due to entrainment is again of a modest scale (Spruit 2015).

5.5 Caveats

The estimates for the long-term effect of entrainment on the growth of convective regions in Section 5.4 are to be taken with caution, however. They are not only crude, time-integrated zeroth-order estimates; the entrainment law (36) is by no means set in stone. Current astrophysical 3D simulations only probe a limited range in the critical parameter Ri_B, and tend to suffer from insufficient resolution for high Ri_B, as shear instabilities develop on smaller and smaller scales.

As a result, it cannot be excluded that the entrainment law (35) transitions to a steeper slope in the astrophysically relevant regime of high Ri_B. Experiments also compete with the difficulties of a limited dynamic range in Reynolds, Prandtl, and Péclet number, and remain inconclusive about the regime of high Ri_B that obtains in stellar interiors. Power-law exponents larger than n = 1 (up to n = 7/4) have also been reported in this regime as alternatives to n = 1 (Fernando 1991; Strang & Fernando 2001; Fedorovich, Conzemius, & Mironov 2004). A power-law exponent n > 1 would imply a strong suppression of entrainment in stellar interiors under most circumstances, and the long-term effect of entrainment would be negligible. Moreover, magnetic fields will affect the shear-driven instabilities responsible for convective boundary mixing (Brüggen & Hillebrandt 2001).

Finally, most of the current 3D simulations of convective boundary mixing suffer from another potential problem; the balance between nuclear energy generation and neutrino cooling that obtains during quasi-stationary shell burning stages is typically violated, or neutrino cooling is not modelled at all. Jones et al. (2016b) pointed out that this may be problematic if neutrino cooling decelerates the buoyant convective plumes and reduces the shear velocity at the interfacial boundary. Only sufficiently long simulations will be able clarify whether the strong entrainment seen in some numerical simulations is robust or (partly) specific to a transient adjustment phase.

Thus, it remains to be seen whether convective boundary mixing has significant effects on the structure of supernova progenitors. Even if it does, it is not clear whether it will qualitatively affect the landscape of supernova progenitors. The general picture of the evolution of massive stars may stay well within the bounds of the variations that have been explored already, albeit in a more parametric way (see, e.g., Sukhbold & Woosley 2014).

6 CONCLUSIONS

It is evident that our understanding of the supernova explosion mechanism has progressed considerably over the last few years. Whilst simulations of CCSNe have yet to demonstrate that they can correctly reproduce and explain the whole range explosions that is observed in nature, there are plenty of ideas for solving the remaining problems. Some important milestones from the last few years have been discussed in this paper, and can be summarised as follows:

• • ECSN-like explosions of supernova progenitors with the lowest masses ( $8\ldots 10\,\text{M}_\odot$ ) can be modelled successfully both in 2D and in 3D. Regardless of the precise evolutionary channel from which they originate, supernovae from the transition region between the super-AGB star channel and classical iron-CCSNe share similar characteristics, i.e., low explosion energies of $\mathord {\sim }10^{50} \, {\rm erg}$ and small nickel masses of a few $10^{-3}\,\text{M}_\odot$ . Due to the ejection of slightly neutron-rich material in the early ejecta, they are an interesting source site for the production of the lighter neutron-rich trans-iron elements (Sr, Y, Zr), and are potentially even a site for a weak r-process up to Ag and Pd (Wanajo et al. 2011). An unambiguous identification of ECSN-like explosions amongst observed transients is still pending; however, although there are various candidate events.

• • Though it has yet to be demonstrated that the neutrino-driven explosion mechanism can robustly account for the explosions of more massive progenitors, first successful 3D models employing multi-group neutrino transport have recently become available. The reluctance of the first 3D models to develop explosions due to the different nature of turbulence in 3D proves to be no insurmountable setback; and even the unsuccessful 3D models computed so far appear to be close to explosion.

• • Some of the recent 2D models produced by different groups (Summa et al. 2016; O’Connor & Couch 2015) show similar results, which inspires some confidence that the simulations are now at a stage where modelling uncertainties due to different numerical methodologies are under reasonable control, though they have not been completely eliminated yet. We have addressed some of the sensitivities to the modelling assumption in this paper, including possible effects of numerical resolution as well as various aspects of the neutrino transport treatment.

• • Recent studies have helped to unravel how the interplay between neutrino heating and hydrodynamic instabilities works quantitatively, and they have clarified why neutrino-driven mechanism can be obtained with a considerably smaller driving luminosity in multi-D.

• • There is a number of ideas about missing physics that could make the neutrino-driven mechanism robust for a wider range of progenitors. These include rapid rotation (Nakamura et al. 2014; Janka et al. 2016; though stellar evolution makes this unlikely as a generic explanation), changes in the neutrino opacities (Melson et al. 2015b), and a stronger forcing of non-radial instabilities due to seed perturbations from convective shell burning (Couch & Ott 2013; Couch et al. 2015; Müller & Janka 2015; Müller et al. 2016a).

• • 3D initial conditions for supernova simulations have now become available (Couch et al. 2015; Müller et al. 2016a), and promise to play a significant and beneficial role in the explosion mechanism. A first 3D multi-group simulation starting from a 3D initial model of an $18\,\text{M}_\odot$ progenitor has been presented in this review. The model has already reached an explosion energy of 5 × 10⁵⁰ erg, and suggests that the observed range of explosion energies may be within reach of 3D simulations.

• • Nonetheless, the study of 3D effects in supernova progenitors is yet in its infancy. A thorough exploration of the parameter space is required in order to judge whether they are generically important for our understanding of supernova explosions. This is not only true with regard to the 3D pre-collapse perturbations from shell burning that are crucial to the ‘perturbation-aided’ neutrino-driven mechanism. The role of convective boundary mixing on the structure of supernova progenitors also deserves to be explored.

Many of these developments are encouraging, though there are also hints of new uncertainties that may plague supernova theory in the future. Whether the new ideas of recent years will prove sufficient to explain shock revival in CCSNe remains to be seen. The perspectives are certainly good, but obviously a lot more remains to be done before simulations and theory can fully explain the diversity of core-collapse events in nature. There is no need to fear a shortage of fruitful scientific problems concerning the explosions of massive stars.