Electron-Capture and Low-Mass Iron-Core-Collapse Supernovae: New Neutrino-Radiation-Hydrodynamics Simulations

TL;DR

This study uses new 1D and 2D simulations to explore electron-capture and low-mass iron-core-collapse supernovae, revealing conditions for explosion and the significant role of PNS convection in neutrino emission.

Contribution

The paper provides detailed simulations of various progenitors, demonstrating the ease of explosion in ECSN-like models and highlighting the underestimated impact of PNS convection.

Findings

ECSN and ECSN-like progenitors explode easily in 1D with low energies.

9-11 solar mass progenitors do not explode in 1D and are not easier in 2D.

PNS convection can double neutrino luminosity, affecting supernova dynamics.

Abstract

We present new 1D (spherical) and 2D (axisymmetric) simulations of electron-capture (EC) and low-mass iron-core-collapse supernovae (SN). We consider six progenitor models: the ECSN progenitor from Nomoto (1984, 1987); two ECSN-like low-mass low-metallicity iron core progenitors from Heger (private communication); and the 9-, 10-, and 11-$M_\odot$ (zero-age main sequence) progenitors from Sukhbold et al. (2016). We confirm that the ECSN and ESCN-like progenitors explode easily even in 1D with explosion energies of up to a 0.15 Bethes ($1 {\rm B} \equiv 10^{51}\ {\rm erg}$), and are a viable mechanism for the production of very low-mass neutron stars. However, the 9-, 10-, and 11-$M_\odot$ progenitors do not explode in 1D and are not even necessarily easier to explode than higher-mass progenitor stars in 2D. We study the effect of perturbations and of changes to the microphysics and we find that relatively small changes can result in qualitatively different outcomes, even in 1D, for models sufficiently close to the explosion threshold. Finally, we revisit the impact of convection below the protoneutron star (PNS) surface. We analyze, 1D and 2D evolutions of PNSs subject to the same boundary conditions. We find that the impact of PNS convection has been underestimated in previous studies and could result in an increase of the neutrino luminosity by up to factors of two.

Figures

Tables

Table 1. Table 1 Details of the setup for the perturbed models.

Prog.	$r_1^{-}$	$r_1^{+}$	${\ell }_1$	$n_1$	$\delta v_1$	$r_2^{-}$	$r_2^{+}$	${\ell }_2$	$n_2$	$\delta v_2$	$r_3^{-}$	$r_3^{+}$	${\ell }_3$	$n_3$	$\delta v_3$
	[km]	[km]			$\left[{10}^8\mathrm{cm}{s}^{-1}\right]$	[km]	[km]			$\left[{10}^8\mathrm{cm}{s}^{-1}\right]$	[km]	[km]			$\left[{10}^8\mathrm{cm}{s}^{-1}\right]$
9.0	850	1,350	12	1	1.0	2,500	9,750	10	1	0.3	15,000	450,000	4	1	0.3
10.0	480	0,560	12	1	1.0	1,150	1,800	10	3	1.0	02,500	003,400	4	1	0.5
11.0	450	0,520	12	1	0.5	1,100	1,400	10	1	1.2	01,550	011,000	4	1	0.5

Table 2. Table 2 Summary of models and energy budget in the region V : 100 km ≤ r ≤ 20 , 000 km V:100\ {\rm km}\leq r\leq 20,\!000\ {\rm km} . Values are given at final simulation time.

Prog.	Setup	${10}^3{\xi }_{2.5}$aafootnotemark:	$E_{\mathrm{bind}}$bbfootnotemark:	$U$ccfootnotemark:	$E_0$ddfootnotemark:	$K_r$eefootnotemark:	$K_{\theta }$fffootnotemark:	$E_g$ggfootnotemark:	$E_{\mathrm{tot}}$hhfootnotemark:	${\dot{E}}_{\mathrm{tot}}$iifootnotemark:	$E_{{\nu }_e}$jjfootnotemark:	$E_{{\overline{\nu }}_e}$kkfootnotemark:	$E_{{\nu }_{\mu }}$llfootnotemark:	$t_{\mathrm{end}}$mmfootnotemark:
			$\left[{10}^{50}\mathrm{erg}\right]$	$\left[{10}^{50}\mathrm{erg}\right]$	$\left[{10}^{50}\mathrm{erg}\right]$	$\left[{10}^{50}\mathrm{erg}\right]$	$\left[{10}^{50}\mathrm{erg}\right]$	$\left[{10}^{50}\mathrm{erg}\right]$	$\left[{10}^{50}\mathrm{erg}\right]$	$\left[{10}^{50}\mathrm{erg}{s}^{-1}\right]$	$\left[{10}^{52}\mathrm{erg}\right]$	$\left[{10}^{52}\mathrm{erg}\right]$	$\left[{10}^{52}\mathrm{erg}\right]$	$\left[s\right]$
n8.8	Baseline 1D	-	$0.000$	$0.328$	$-0.058$	$0.994$	$0.000$	$-0.061$	$1.202$	$0.038$	$1.346$	$0.830$	$2.582$	$0.636$
n8.8	ManyBody 1D	-	$0.000$	$0.354$	$-0.062$	$1.158$	$0.000$	$-0.064$	$1.386$	$0.058$	$1.383$	$0.871$	$2.977$	$0.622$
n8.8	Baseline 2D	-	$0.000$	$0.464$	$-0.094$	$1.381$	$0.000$	$-0.105$	$1.646$	$0.265$	$1.555$	$0.918$	$3.150$	$0.619$
n8.8	ManyBody 2D	-	$0.000$	$0.487$	$-0.097$	$1.500$	$0.000$	$-0.107$	$1.783$	$0.287$	$1.566$	$0.930$	$3.412$	$0.602$
u8.1	Baseline 1D	$0.094$	$-0.018$	$0.136$	$-0.076$	$0.016$	$0.000$	$-0.098$	$-0.021$	$0.013$	$1.598$	$1.228$	$3.414$	$0.888$
u8.1	ManyBody 1D	$0.094$	$-0.018$	$0.355$	$-0.115$	$0.063$	$0.000$	$-0.117$	$0.185$	$0.036$	$1.642$	$1.278$	$4.095$	$0.914$
u8.1	Baseline 2D	$0.094$	$-0.018$	$0.645$	$-0.204$	$0.523$	$0.001$	$-0.105$	$0.861$	$0.085$	$2.037$	$1.526$	$5.145$	$1.095$
u8.1	ManyBody 2D	$0.094$	$-0.018$	$0.715$	$-0.211$	$0.636$	$0.001$	$-0.112$	$1.029$	$0.113$	$2.039$	$1.523$	$5.557$	$1.022$
z9.6	Baseline 1D	$0.076$	$-0.008$	$0.173$	$-0.052$	$0.076$	$0.000$	$-0.055$	$0.143$	$0.019$	$1.518$	$1.157$	$3.374$	$0.900$
z9.6	ManyBody 1D	$0.076$	$-0.008$	$0.244$	$-0.073$	$0.256$	$0.000$	$-0.046$	$0.381$	$0.021$	$1.593$	$1.235$	$4.072$	$0.916$
z9.6	Baseline 2D	$0.076$	$-0.008$	$0.472$	$-0.128$	$0.785$	$0.001$	$-0.091$	$1.038$	$0.147$	$1.737$	$1.231$	$4.070$	$0.771$
z9.6	ManyBody 2D	$0.076$	$-0.008$	$0.527$	$-0.137$	$0.931$	$0.001$	$-0.101$	$1.221$	$0.218$	$1.726$	$1.220$	$4.303$	$0.722$
9.0	Baseline 1D	$0.038$	$-0.021$	$0.401$	$-0.206$	$0.034$	$0.000$	$-0.269$	$-0.040$	$0.012$	$1.652$	$1.298$	$3.287$	$0.816$
9.0	ManyBody 1D	$0.038$	$-0.021$	$0.438$	$-0.208$	$0.027$	$0.000$	$-0.261$	$-0.003$	$0.003$	$1.744$	$1.392$	$3.974$	$0.848$
9.0	Baseline 2D	$0.038$	$-0.021$	$0.717$	$-0.273$	$0.168$	$0.014$	$-0.218$	$0.409$	$0.094$	$2.326$	$1.822$	$5.944$	$1.454$
9.0	Baseline Perturb 2D	$0.038$	$-0.021$	$0.728$	$-0.289$	$0.232$	$0.013$	$-0.189$	$0.495$	$0.070$	$2.427$	$1.925$	$6.371$	$1.645$
9.0	ManyBody 2D	$0.038$	$-0.021$	$0.841$	$-0.292$	$0.295$	$0.012$	$-0.229$	$0.627$	$0.202$	$2.361$	$1.846$	$6.655$	$1.381$
10.0	Baseline 1D	$0.216$	$-0.095$	$1.151$	$-0.578$	$0.225$	$0.000$	$-0.996$	$-0.197$	$0.563$	$2.279$	$1.842$	$3.744$	$0.807$
10.0	ManyBody 1D	$0.216$	$-0.095$	$1.226$	$-0.560$	$0.184$	$0.000$	$-0.904$	$-0.053$	$0.477$	$2.458$	$2.016$	$4.627$	$0.875$
10.0	Baseline 2D	$0.216$	$-0.095$	$0.788$	$-0.394$	$0.054$	$0.000$	$-0.370$	$0.078$	$-0.012$	$3.764$	$3.142$	$8.968$	$2.147$
10.0	Baseline Perturb 2D	$0.216$	$-0.095$	$0.996$	$-0.463$	$0.086$	$0.000$	$-0.592$	$0.027$	$0.186$	$3.059$	$2.450$	$6.505$	$1.272$
10.0	ManyBody Perturb 2D	$0.216$	$-0.095$	$1.628$	$-0.593$	$0.293$	$0.024$	$-0.591$	$0.762$	$0.425$	$3.245$	$2.645$	$8.413$	$1.575$
10.0	ManyBody 2D	$0.216$	$-0.095$	$1.126$	$-0.482$	$0.087$	$0.000$	$-0.615$	$0.117$	$0.303$	$3.082$	$2.463$	$7.183$	$1.216$
11.0	Baseline 1D	$7.669$	$-0.170$	$2.794$	$-1.530$	$0.363$	$0.000$	$-2.188$	$-0.560$	$0.692$	$2.449$	$2.002$	$3.942$	$0.924$
11.0	ManyBody 1D	$7.669$	$-0.170$	$3.019$	$-1.536$	$0.339$	$0.000$	$-2.174$	$-0.351$	$0.781$	$2.558$	$2.105$	$4.724$	$0.927$
11.0	Baseline 2D	$7.669$	$-0.170$	$3.751$	$-1.573$	$0.700$	$0.068$	$-1.628$	$1.318$	$0.497$	$3.136$	$2.539$	$7.033$	$1.508$
11.0	Baseline Perturb 2D	$7.669$	$-0.170$	$3.585$	$-1.544$	$0.379$	$0.068$	$-1.816$	$0.671$	$0.619$	$3.311$	$2.707$	$7.428$	$1.647$
11.0	ManyBody 2D	$7.669$	$-0.170$	$3.929$	$-1.674$	$0.578$	$0.040$	$-1.521$	$1.352$	$1.089$	$3.370$	$2.761$	$8.986$	$1.780$

Table 3. Table 3 PNS star masses at the final simulation time.

Prog.	Setup	$M_{\mathrm{Baryon}}$aafootnotemark:	$M_{\mathrm{Grav}}$bbfootnotemark:	${\dot{M}}_{\mathrm{Baryon}}$ccfootnotemark:
		$\left[M_{\odot }\right]$	$\left[M_{\odot }\right]$	$\left[M_{\odot }{s}^{-1}\right]$
n8.8	Baseline 1D	$1.300$	$1.193$	$0.003$
n8.8	ManyBody 1D	$1.299$	$1.193$	$0.003$
n8.8	Baseline 2D	$1.294$	$1.188$	$0.004$
n8.8	ManyBody 2D	$1.294$	$1.188$	$0.005$
u8.1	Baseline 1D	$1.391$	$1.270$	$0.005$
u8.1	ManyBody 1D	$1.385$	$1.265$	$0.000$
u8.1	Baseline 2D	$1.374$	$1.256$	$0.008$
u8.1	ManyBody 2D	$1.372$	$1.254$	$0.000$
z9.6	Baseline 1D	$1.379$	$1.260$	$0.000$
z9.6	ManyBody 1D	$1.376$	$1.257$	$0.000$
z9.6	Baseline 2D	$1.367$	$1.250$	$0.002$
z9.6	ManyBody 2D	$1.366$	$1.249$	$0.003$
9.0	Baseline 1D	$1.369$	$1.252$	$0.013$
9.0	ManyBody 1D	$1.369$	$1.252$	$0.011$
9.0	Baseline 2D	$1.361$	$1.245$	$0.000$
9.0	Baseline Perturb 2D	$1.359$	$1.243$	$0.001$
9.0	ManyBody 2D	$1.359$	$1.243$	$0.001$
10.0	Baseline 1D	$1.509$	$1.368$	$0.067$
10.0	ManyBody 1D	$1.512$	$1.371$	$0.056$
10.0	Baseline 2D	$1.540$	$1.394$	$0.008$
10.0	Baseline Perturb 2D	$1.528$	$1.385$	$0.024$
10.0	ManyBody Perturb 2D	$1.510$	$1.369$	$0.003$
10.0	ManyBody 2D	$1.525$	$1.382$	$0.026$
11.0	Baseline 1D	$1.515$	$1.373$	$0.093$
11.0	ManyBody 1D	$1.513$	$1.372$	$0.089$
11.0	Baseline 2D	$1.514$	$1.372$	$0.034$
11.0	Baseline Perturb 2D	$1.519$	$1.377$	$0.024$
11.0	ManyBody 2D	$1.496$	$1.358$	$0.024$

Equations

$$\delta\mathbf{v}_{i}=\left\{\begin{array}[]{ll}\frac{C_{i}}{\rho}\nabla\times\bm{\Psi}_{i},&r_{i}^{-}\leq r\leq r_{i}^{+},\\ 0,&\textrm{otherwise};\end{array}\right.$$

$${\mathbf{\Psi }}_i={\mathbf{e}}_{\varphi }\frac{\sqrt{\sin \theta }}{r}\sin \left(n_i\pi \frac{r-r_i^{-}}{r_i^{+}-r_i^{-}}\right)Y_{{\ell }_i,1}(\theta ,0)$$

Full text

††slugcomment: Draft version

1. Introduction

The formation of a massive iron core at the end of the evolution of stars with zero-age main-sequence (ZAMS) masses larger than ${\sim}12\ M_{\odot}$ is a robust prediction of stellar evolution theory. These stars undergo core-collapse once their cores reach the Chandrasekhar mass and may explode as core-collapse supernovae.

The fate of less massive stars in the ZAMS range ${\sim}8\ M_{\odot}$ to ${\sim}12\ M_{\odot}$ is less clear. Depending on the initial mass, the ultimate fate could be to form massive white dwarfs, to form iron cores as do regular massive stars, or to explode as electron-capture supernovae before forming an iron core (e.g. Nomoto, 1984, 1987; Jones et al., 2013; Doherty et al., 2015; Woosley & Heger, 2015; Doherty et al., 2017). It has also been suggested that these stars might undergo violent flashes and power unusual transients before their deaths (Woosley & Heger, 2015; Jones et al., 2016).

ECSNe and low-mass iron-core CCSNe with similar features are expected to occur in a relatively narrow range of ZAMS masses. However, they might account for a significant fraction of gravitational-collapse supernovae, given that the initial mass function of stars drops rapidly towards high masses. These progenitors have compact cores with tenuous envelopes, which result in a steep drop of the accretion rate after core bounce. This, in turn, triggers early explosions, even under the assumption of spherical symmetry (Kitaura et al., 2006; Janka et al., 2008; Burrows et al., 2007; Fischer et al., 2010; Janka et al., 2012). ECSNe and ECSN-like CCSNe are expected to be underenergetic and possibly underluminous and to have small ${}^{56}{\rm Ni}$ yields and peculiar nucleosynthetic abundances (e.g. Nomoto et al., 1982; Kitaura et al., 2006; Janka et al., 2008; Hoffman et al., 2008; Wanajo et al., 2011; Melson et al., 2015b; Wanajo et al., 2017).

Kitaura et al. (2006) suggested that ECSN-like events might explain a subclass of Type-II SNe with unusually low luminosities (Pastorello et al., 2004; Spiro et al., 2014). An ECSN has also been invoked to explain SN 1054 and the associated Crab remnant (Nomoto et al., 1982; Takahashi et al., 2013; Smith, 2013; Tominaga et al., 2013). According to historical records, SN 1054 was not underluminous. However, SN-1054 was likely underenergetic, with an explosion energy around $10^{50}$ erg, as indicated by the low-mass of the Crab nebula’s filaments and their relatively low expansion velocity, as well as by the small inferred ${}^{56}{\rm Ni}$ yield (Müller, 2016, and references therein). On the basis of measured isotopic abundance anomalies, it has also been suggested that a low-mass CCSN might have been the trigger that started the formation of our Solar System (Banerjee et al., 2016).

ECSNe and ECSN-like progenitors have attracted significant interest in the CCSN mechanism community due to their “explodability” and the fact that they allow for self-consistent studies also in 1D (with the assumption of spherical symmetry). Hillebrandt et al. (1984) performed the first 1D simulations of the collapse, bounce, and explosion of the original n8.8 ECSN progenitor of Nomoto (Nomoto, 1984, 1987), using an approximate gray neutrino transport scheme. They found an energetic explosion by the prompt shock mechanism with an energy of ${\sim}2\cdot 10^{51}$ erg. Subsequent studies, with modern neutrino interactions and multi-group transport in 1D and 2D, performed by Kitaura et al. (2006); Janka et al. (2008); Burrows et al. (2007) and Fischer et al. (2010) found much weaker explosions (${\sim}10^{50}\ {\rm erg}$) powered by the delayed neutrino mechanism. Müller et al. (2012a) considered an $8.1\ M_{\odot}$ (ZAMS) progenitor with metallicity $Z=10^{-4}$, u8.1, which formed an iron core, but had a stellar structure very similar to the n8.8 progenitor, and found a similarly early explosion. Another iron-core progenitor, the zero-metallicity $9.6$-$M_{\odot}$ model, z9.6, was considered and determined to have a qualitatively similar outcome by Janka et al. (2012); Müller et al. (2013); Müller & Janka (2014) in 2D, and by Melson et al. (2015b) in 3D. Very recently, Wanajo et al. (2017) presented new nucleosynthetic calculations for the n8.8, u8.1, and z9.6 models, as well as a summary of their associated explosion characteristics with the Coconut-Vertex code.

The low-mass, but otherwise “canonical,” $11.2$-$M_{\odot}$ progenitor from Woosley et al. (2002) and the $12.0$-$M_{\odot}$ progenitor from Woosley & Heger (2007) have been considered by several groups, (e.g. Buras et al., 2006; Takiwaki et al., 2012; Müller et al., 2012b; Bruenn et al., 2013, 2016; Dolence et al., 2015; Müller, 2015; Summa et al., 2016; O’Connor & Couch, 2015; Burrows et al., 2016; Nagakura et al., 2017). While the $11.2$-$M_{\odot}$ progenitor has a post-bounce evolution that is qualitatively similar to the u8.1 progenitor (Müller et al., 2012a), the $12.0$-$M_{\odot}$ progenitor either explodes very late (Summa et al., 2016) or not at all (O’Connor & Couch, 2015; Dolence et al., 2015; Skinner et al., 2016), with the exception of the simulations by Bruenn et al. (2013, 2016). This challenges the commonly held idea that low-mass iron-core progenitors should explode easily and similarly to an ECSN.

With the goal of characterizing the different explosions of ECSNe, ECSN-like CCSN, and canonical, but low-mass, CCSNe, we present 1D and 2D Fornax simulations (Dolence et al. 2017, in prep.) of the collapse, bounce, and subsequent evolution of six progenitor models. We consider the ECSN n8.8 progenitor from Nomoto (1984, 1987), the ECSN-like u8.1 and z9.6 progenitors from Heger (private communication), and the 9-, 10-, and $11$-$M_{\odot}$ solar-metallicity iron-core collapse progenitors from Sukhbold et al. (2016). We show that low-mass CCSNe always fail to explode in 1D and evolve in a qualitatively different way compared to ECSNe and ECSN-like CCSNe. We also discuss the impact of the dimensionality, of pre-supernova perturbations, and of changes in the microphysics. In particular, for the latter, we focus on the effects of many-body corrections to the axial-vector term in the neutrino-nucleon scattering rate recently explored by Horowitz et al. (2017).

The rest of this paper is organized as follows. First, in Section 2, we give an overview of the simulation setup and of the properties of the progenitor models. We discuss the qualitative outcome of our simulations in Section 3, while a more quantitative account of the energetics of the explosions is given in Section 4. We discuss the properties of the neutrino radiation in Section 5. Section 6 is dedicated to the properties and evolution of the remnant protoneutron stars. Finally, we summarize and discuss our results in Section 7.

2. Progenitors and Setup

As previously discussed, we consider six progenitor models, which we label as n8.8 (Nomoto, 1984, 1987), u8.1 and z9.6 (Heger private communication), and 9.0-, 10.0-, $11.0$-$M_{\odot}$ (Sukhbold et al., 2016). The u8.1 progenitor has metallicity $10^{-4}$ of Solar, the z9.6 has zero metallicity, and all other progenitors have Solar metallicity. All of the progenitors have been evolved up to the point of core-collapse, defined as the time their radial infall velocity has reached ${\sim}1000$ ${\rm km}\ {\rm s}^{-1}$. Their structure (density and exterior binding energies) are shown in Fig. 1. All of these progenitors have relatively compact cores and loosely bound envelopes. The values of the compactness parameter $\xi_{2.5}$ (O’Connor & Ott, 2011, 2013) computed from the progenitor models are given in Tab. 2. They range from $\simeq 7.6\cdot 10^{-5}$, for the z9.6 progenitor, to $\simeq 7.7\cdot 10^{-3}$, for the 11.0-$M_{\odot}$ progenitor. $\xi_{2.5}$ cannot be computed for the n8.8 progenitor since it only extends to $\simeq 1.32\ M_{\odot}$. Note that, for numerical reasons, we modify the n8.8 progenitor for $\rho\leq 10^{4}\ {\rm g}\cdot{\rm cm}^{-3}$ with the addition of a constant temperature envelope with $\rho\propto r^{-2}$. We verified, by changing the density of the envelope over 3 orders of magnitude, that the explosion energy of the n8.8 progenitor is not very sensitive to this envelope, although the shock propagation speed is obviously affected.

All of the progenitors are non-rotating and have been evolved in 1D. There is currently significant interest in the possible impact of pre-supernova turbulence on the development of explosions, which several authors have determined to be beneficial and, in some cases, crucial to the outcome (Couch & Ott, 2013; Müller & Janka, 2015; Abdikamalov et al., 2016; Takahashi et al., 2016; Müller et al., 2017). Of relevance, the first progenitor models evolved in 3D shortly before core collapse have also recently become available (Couch et al., 2015; Müller et al., 2016). In this spirit, we also consider the impact of perturbations on the 9.0-, 10.0-, and 11.0-$M_{\odot}$ progenitors using the approach introduced by Müller & Janka (2015), which we briefly describe. We introduce velocity perturbations in three regions $r_{i}^{-}\leq r\leq r_{i}^{+}$, $i=1,2,3$, obeying the divergence-free condition $\nabla\cdot(\rho\,\delta\mathbf{v}_{i})=0$. These are generated as

[TABLE]

where

[TABLE]

and $n_{i}$, $\ell_{i}$ are the number of convective cells in the radial and angular directions, respectively. Finally, $C_{i}$ is tuned to achieve a given maximum perturbation amplitude. The parameters we use are given in Tab. 1.

We evolve these progenitors from the onset of collapse with the neutrino-radiation-hydrodynamics code Fornax (Skinner et al. 2016; Burrows et al. 2016; Dolence et al. 2017, in prep.). Fornax solves for the transport of neutrinos using a multi-dimensional moment scheme with an analytic closure for the \nth2 and \nth3 moments (Shibata et al., 2011; Murchikova et al., 2017). Similar moment methods have also been recently adopted for the CCSN problem by other groups (e.g. O’Connor, 2015; Just et al., 2015; O’Connor & Couch, 2015; Roberts et al., 2016). The moment equations are solved using a \nth2-order finite-volume scheme with the HLLE approximate Riemann solver (Einfeldt, 1988), modified as in Audit et al. (2002) and O’Connor (2015) to reduce the numerical dissipation in the diffusive limit. Fornax separately evolves electron neutrinos $\nu_{e}$ and anti-electron neutrinos $\bar{\nu}_{e}$, while heavy-lepton neutrinos $\nu_{\mu}$, $\nu_{\tau}$, and the respective anti-particles are lumped together as a single species, which we denote as “$\nu_{\mu}$” (see Bollig et al. 2017 for a discussion of the possible limitation of this aporach). The energy spectra of neutrinos are resolved using 20 logarithmically-spaced energy groups extending to $300\ {\rm MeV}$ for electron neutrinos and to $100\ {\rm MeV}$ for anti-electron and heavy-lepton neutrinos.

The set of neutrino-matter interactions included in our simulations are described Burrows et al. (2006). We include weak magnetism and recoil correction to neutrino-nucleon scattering and absorption (Horowitz, 2002). We treat inelastic neutrino-electron scattering with the scheme of Thompson et al. (2003) and the relativistic formalism summarized in Reddy et al. (1999). Inelastic neutrino-nucleon scattering is included using the formalism of Thompson et al. (2003), for $\nu_{e}$ and $\bar{\nu}_{e}$, while we use the approach of Müller & Janka (2015) for inelastic scattering of heavy-lepton neutrinos on nucleons. For the latter, we use $6k_{\rm B}T$, instead of $3k_{\rm B}T$ for the crossover energy between upscattering and downscattering (Thompson et al., 2000; Tubbs, 1979), where $k_{\rm B}$ is Boltzmann’s constant and $T$ is the temperature. That is, we approximate the redistribution rate of the heavy-lepton neutrinos to be proportional to $(\epsilon_{\nu}-6k_{\rm B}T)/m_{n}c^{2}$, where $m_{n}$ is the neutron mass and $c$ is the speed of light, and $\epsilon_{\nu}$ is the incoming neutrino energy. Electron capture on heavy nuclei during the infall is treated following Bruenn (1985). However, we disable electron capture on heavy nuclei for the n8.8 progenitor to prevent an unphysical neutronization burst during the infall, which is due to the failure of our nuclear statistical equilibrium approximation in the outer core of this progenitor.

We perform two variants of each simulation. “Baseline” includes all the neutrino-matter interaction discussed above. “ManyBody” also include many-body corrections to the neutrino-nucleon scattering cross-section as estimated by Horowitz et al. (2017). These are implemented using the fit to the axial response factor $S_{A}$ they provided. The runs with perturbations are performed using the Baseline physics setup, with the exception the 10.0-$M_{\odot}$ which is evolved with both the Baseline and the ManyBody setup.

The hydrodynamic equations are solved using a high-resolution shock-capturing scheme with \nth3-order reconstruction and the HLLC approximate Riemann solver (Toro et al., 1994). The details of the numerical schemes are discussed in (Skinner et al. 2016; Burrows et al. 2016; Dolence et al. 2017, in prep.). For the simulations presented here, we use a spherical grid with 678 points extending up to 20,000 km. The grid has a constant spacing $\Delta r$ of 0.5 km for $r\lesssim 10\ {\rm km}$ and then smoothly transitions to a logarithmically spaced grid with $\Delta r/r\simeq 0.01$ for $r\gtrsim 100\ {\rm km}$. For the 2D simulations, we use 256 angular zones with angular resolution smoothly varying between $\simeq 0.95^{\circ}$ at the poles and $\simeq 0.64^{\circ}$ at the equator. The angular grid is also progressively derefined towards the center, *i.e.*we use a dendritic grid (Dolence et al. 2017, in prep.), to avoid an excessively restrictive CFL condition in the angular direction.

We adopt the Lattimer-Swesty equation of state with nuclear compressibility parameter $220\ {\rm MeV}$ (Lattimer & Swesty, 1991) and treat gravity in the monopole approximation using a general-relativistic potential, following Marek et al. (2006).

Finally, we carry out all 2D simulations until the maximum shock radius exceeds 19,000 km, or until the explosion is deemed unsuccessful.

3. Overall Dynamics

A first glance of our results can be gained from Fig. 2, which shows the average shock radii for all progenitors in 1D and 2D with both the Baseline and ManyBody setups. As in previous works by others (Kitaura et al., 2006; Janka et al., 2008; Burrows et al., 2007; Fischer et al., 2010; Müller et al., 2012a; Janka et al., 2012; Müller et al., 2013; Melson et al., 2015b), we find early explosions for the n8.8, z9.6, and u8.1 progenitors. The same progenitors also explode in 1D spherical symmetry, although the u8.1 only with the ManyBody setup.

In the n8.8 models, the shock suddenly accelerates outwards starting from ${\sim}50\ {\rm ms}$ after bounce, a few tens of milliseconds earlier than in Kitaura et al. (2006) and Janka et al. (2008), but similarly to Fischer et al. (2010). We note that these works used different equations of state compared to us. Kitaura et al. (2006) and Janka et al. (2008) used the Lattimer-Swesty equation of state with nuclear compressibility parameter $180\ {\rm MeV}$ (Lattimer & Swesty, 1991), while Fischer et al. (2010) used the equation of state of Shen et al. (1998). These studies also differ in their treatment of nuclear burning and electron capture during infall. These differences result in variations of the positions at which the shock originally forms, which in turn might explain the different explosion times. In the simulations of the Garching group the homologous core at bounce has a mass of $0.425\ M_{\odot}$ (Kitaura et al., 2006), while Fischer et al. (2010) report a larger core mass of $0.625\ M_{\odot}$. In our simulations, the homologous core encloses a mass of $0.611\ M_{\odot}$ at bounce.

None of the progenitors from Sukhbold et al. (2016) explode in self-consistent 1D simulations. Somewhat surprisingly, the 10.0-$M_{\odot}$ progenitor fails to explode, within the simulation time, also in 2D. However, both the 9.0- and 11.0-$M_{\odot}$ progenitor explode successfully in 2D, with either the Baseline or the ManyBody setups.

To visualize the different outcomes of the progenitors from Sukhbold et al. (2016) in 2D, we highlight their early post-bounce average shock radii in Fig. 3. The 9.0-$M_{\odot}$ and 11.0-$M_{\odot}$ models have delayed explosions at ${\sim}0.3\ {\rm s}$ after bounce (Baseline setup) or ${\sim}0.2\ {\rm s}$ after bounce (11.0-$M_{\odot}$, with ManyBody). The inclusion of perturbations does not affect the outcome of the 9.0-$M_{\odot}$ progenitor. Surprisingly, perturbations result in a somewhat weaker explosion for the 11.0-$M_{\odot}$ progenitor, as can be inferred from the smaller shock expansion velocity and as quantified in Sec. 4. In the case of the 10.0-$M_{\odot}$ progenitor, perturbations are able to trigger a weak explosion ${\sim}0.4\ {\rm s}$ after bounce, but only in combination with the many-body effects included in the ManyBody setup. This is the only model for which we find perturbations to yield a qualitative change to the evolution, despite the fact that the amplitude of the initial perturbations is at the upper end of what could be considered as realistic, with turbulent velocities reaching ${\sim}1000\ {\rm km}\ {\rm s}^{-1}$.

Some insight into the reason for the different evolutions can be gained from the analysis of the accretion rate history of the progenitors, as recently suggested by Suwa et al. (2016) and Müller (2016). Fig. 4 shows the accretion rate at 500 km for all progenitors in 1D, with the Baseline setup. Since the 1D progenitors from Sukhbold et al. (2016) fail to explode, these are “intrinsic” accretion rates, not affected by the explosion. The accretion rates for the n8.8 and z9.6 are unaffected by the explosions up to the point where they are shown (afterwards, they become negative as the inflow turns into an outflow). It is easily seen that the accretion rates for progenitors exploding in 1D decline steeply at very early times, which sets them apart from “normal” massive stars, as also pointed out by Müller (2016). We find that the accretion rate of the 10.0-$M_{\odot}$ progenitor is significantly higher than for the 9.0- and 11.0-$M_{\odot}$ models during the critical phase when the other two start exploding. The sudden growth of the accretion rate of the 10.0-$M_{\odot}$ progenitor around ${\sim}0.2\ {\rm s}$ after bounce is due to a small density inversion present in the original progenitor profile (Fig. 1). Obviously, such a density inversion would be Rayleigh-Taylor unstable and is not expected to be present in Nature.

Figures 5 and 6 show snapshots of the entropy and of the pressure contrast, defined following Fernández & Thompson (2009) to be $r\,|\nabla p|/p$, for selected progenitors at $0.15\ {\rm s}$, $0.2\ {\rm s}$, and $0.4\ {\rm s}$ after bounce. In all our models, neutrino-driven convection is seeded in the region behind the shock by perturbations induced by the dendritic grid and develops rapidly following bounce (i.e., starting from ${\sim}0.15\ {\rm s}$ after bounce). Large-scale shock-sloshing motions, possibly due to the standing accretion shock instability (SASI) (Blondin et al., 2003; Foglizzo et al., 2007), are only present at late times in models that fail to explode (e.g., the 10.0-$M_{\odot}$), after the shock has receded to less than $100\ {\rm km}$ in radius. Due to the low post-bounce accretion rates, convection is the dominant instability at early times and for all exploding models (Foglizzo et al., 2006; Burrows et al., 2012; Murphy et al., 2013; Müller et al., 2012a; Ott et al., 2013; Couch & O’Connor, 2014; Fernández et al., 2014; Abdikamalov et al., 2015).

That said, convection appears to have a different role in the onset of the explosion of the n8.8 compared to that of the other progenitors (for the setups resulting in explosions). For the n8.8, the shock starts expanding rapidly already after about one convective overturn, before the high-entropy plumes are able to reach it (Figs. 5 and 6). The onset of this explosion is essentially spherically symmetric and does not seem to be significantly aided by convection, although there are small differences between 1D and 2D visible in Fig. 2.

In the case of the ECSN-like z9.6 and u8.1 progenitors, convective plumes are able to reach the shock and appear to have an important role in triggering the explosion, which is otherwise delayed or fails altogether in 1D. The role of convective instabilities in the explosion of the z9.6 progenitor was also studied by Melson et al. (2015b). Our results are in qualitative agreement with theirs, but there are quantitative differences. In particular, they found a significantly more delayed explosion in 1D than in 2D and 3D, as a comparison of our 1D and 2D shock radius evolutions to theirs (Fig. 3 of Melson et al. 2015b) demonstrates.

Differently from the ECSN and ECSN-like progenitors, the 9.0-, 10.0-, 11.0-$M_{\odot}$ progenitors explode only after a few convective overturns with the emergence of one or more large plumes that succeed in pushing the shock to a sufficiently large radius to trigger a run-away expansion. This behavior is commonly observed in 2D CCSN simulations (e.g., Fernández et al., 2014).

After the explosion sets in, the shock and the material immediately behind it expand almost self-similarly, with roughly constant velocities. Behind them, we observe the emergence of a higher-entropy neutrino-driven wind with dynamics similar to the one reported by Burrows (1987) and Burrows et al. (1995). For the n8.8, z9.6, and u8.1 progenitors the wind is quasi-spherical, it produces weak shocks visible in the pressure-contrast visualizations in Figs. 5 and 6, and drives Rayleigh-Taylor instabilities as it pushes on the slower, heavier material above. In the case of asymmetric explosions, the wind is typically confined in ${\sim}90^{\circ}$ wedges along the axis where it drives the inflation of a large bubble, while fall-back accretion continues along the equator. For these models, late-time accretion is primarily responsible for the growth of the explosion energy, while for the ECSN and ECSN-like explosions the explosion energy injection is due to the wind. We caution the reader, however, that the degree of asymmetry in the 9.0-, 10.0-, and 11.0-$M_{\odot}$ progenitor explosions is likely to be artificially magnified by the assumption of axisymmetry, and we speculate that the neutrino-driven wind will be closer to spherical in full-3D simulations.

{turnpage}

As a representative example of an ECSN-like explosion, we show in Figs. 7 and 8 a summary of the evolution of the n8.8 progenitor in 1D and 2D, evolved with the Baseline setup. Despite their relatively similar average shock trajectories and the fact that the n8.8 shock remains nearly spherical in both 1D and 2D simulations, there are substantial differences between the 1D and 2D explosions. In particular, the entropy in the 1D model exceeds that of the 2D model, while the velocity of the neutrino-driven wind in the 2D model exceeds that in the 1D model, roughly by a factor of two. This is partly due to the tendency of 1D models to create low-density regions that become overheated by neutrinos. More importantly, starting from ${\sim}0.2\ {\rm s}$ after bounce, the 2D simulation shows significantly larger neutrino luminosities (see Sec. 6 for a detailed account). This results in stronger neutrino-driven winds, with increased velocities (see Fig. 8) and correspondingly smaller expansion timescales.

The 9.0-, 10.0-, and 11.0-$M_{\odot}$ progenitors evolve in a qualitatively different way. Fig. 9 shows the evolution in 2D of the 9.0-$M_{\odot}$ progenitor with the Baseline setup, which we take as representative of regular (successful) CCSNe from a low-mass progenitor. The explosion of the 9.0 follows ${\sim}100$ ms of shock stagnation and is triggered at the time when the Si/O interface is accreted. As already mentioned, the explosion is asymmetric. It is also marginal, with small velocities immediately below the shock and a degree of sustained fallback at the end of the simulation. The z9.6 and u8.1 progenitors’ dynamics are similar to that of the regular CCSN progenitors at early times, during the shock stagnation epochs. However, their explosion is also triggered by the sharp accretion rate drop in an almost spherical way, and then evolves in a qualitatively similar way to that of the n8.8 progenitor.

4. Explosion Energetics

We estimate explosion energies using a fixed-volume energy analysis similar to that in Bruenn et al. (2016). We consider the region $r\geq 100\ {\rm km}$, and we compute the integrated total energy $E_{\rm tot}$ as the sum of internal $U$, kinetic $K$, and gravitational binding energy $E_{g}$. To this, we subtract the zero-temperature energy of the material $E_{0}$ computed from the equation of state (EOS) used for the evolution. We remark that $E_{\rm tot}$ includes the contribution of the binding energy of the material exterior to the shock but interior to the computational domain. The values of these quantities at the end of our simulations are given in Tab. 2. There, we also quote the total energy liberated in each species of neutrinos $E_{\nu_{e}}$, $E_{\bar{\nu}_{e}}$, and $E_{\nu_{\mu}}$. The final explosion energy can be estimated by subtracting in absolute value the binding energy of the material exterior to $20,\!000\ {\rm km}$ from $E_{\rm tot}$. This is shown, as a function of time, in Fig. 10. Note that, at late times, the net energy in the integration region can take (small) positive values also for failing models, e.g., for the 10.0-$M_{\odot}$ progenitor. This is mostly because our metric includes with a positive sign also the kinetic energy of infalling fluid elements. Furthermore, we remark that, because of the general relativistic (GR) corrections included in our treatment of gravity, the total energy is not conserved. Instead, the gravitational potential decreases by several percent over the duration of our simulations, because of neutrino losses.

We find explosion energies ranging from a few percent of a Bethe ($1\ {\rm B}\equiv 10^{51}\ {\rm erg}$), like the 9.0-Baseline 2D run, to values in excess of $0.17$ Bethes, for the n8.8-ManyBody 2D simulation. The explosion energy we estimate for the z9.6 progenitor with the Baseline setup is ${\sim}50\%$ larger than that reported by Melson et al. (2015b) for their 2D Prometheus-Vertex simulation. The estimated explosion energies for the z9.6 and u8.1 in 2D are also similarly larger than those of Coconut-Vertex in 2D, as quoted in Wanajo et al. (2017). The discrepancy is somewhat larger for the n8.8 progenitor, where we estimate an explosion energy almost a factor two larger than reported by Wanajo et al. (2017). On the other hand, the explosion energies for the z9.6 and n8.8 in 1D ($0.01\ {\rm B}$ and $0.12\ {\rm B}$) are in good agreement with those reported by the Garching group (Melson et al., 2015b; Kitaura et al., 2006), suggesting that the discrepancies might be due to multi-dimensional effects. Note, however, that in Janka et al. (2008) the Garching group reported a $20\%$ smaller explosion energy for the n8.8 progenitor compared to our results and their own previous calculations (Kitaura et al., 2006). Explosion energies for the progenitors from Sukhbold et al. (2016) have not been reported before, so no comparison is possible.

The reduction of neutral current interactions in the ManyBody setup yields an increase in the explosion energies of 10% to 50%, depending on the model. The amplification is particularly large for the 9.0-$M_{\odot}$ progenitor in 2D and the z9.6 progenitor in 1D, where the ManyBody setup boosts the explosion energy by ${\sim}50\%$. These large amplifications are due to the proximity of these progenitors to criticality, which amplifies their sensitivity to relatively small changes in the input microphysics.

The role of perturbations on the explosion energy (Tab. 2) is not completely clear. In the case of the 9.0- and 10.0-$M_{\odot}$ progenitors, the inclusion of perturbations is beneficial. The first explodes with slightly larger (${\sim}10\%$) explosion energy than without perturbations. The second goes from a failed to a successful, albeit underenergetic, explosion with the introduction of perturbations in combination with many body corrections to neutral current interactions. Somewhat surprisingly, in the case of the 11.0-$M_{\odot}$ progenitor, the inclusion of perturbations results in a reduction of the explosion energy by a factor 2. The reason is that 11.0-Perturb explosion entrains more bound mass and the ejecta lose more energy, while doing work on the infalling envelope of the star. On the other hand, note that the explosion energies for the 11.0-$M_{\odot}$ progenitor are still growing significantly at the end of our simulations, so the difference between the 11.0-Perturb and 11.0-Baseline models might be only transitory.

It is important to keep in mind that the explosion energies we quote are not final. Indeed the estimated explosion energy is still growing significantly at the time we stop the calculation for many of our simulations. This is not surprising in the light of the results of Müller (2015), who studied the development of an explosion in the 11.2-$M_{\odot}$ progenitor from Woosley et al. (2002) and found the explosion energy to saturate only after several seconds. On the other hand, the explosion energies saturate very rapidly for the n8.8, u8.1, and z9.6 progenitors and appear to have converged within the simulation time.

Another caveat is that our estimate of the explosion energy is more conservative than the commonly used “diagnostic explosion energy.” The former is computed as $E_{\rm tot}$, but only integrated over unbound and/or radially expanding fluid elements (e.g., Buras et al., 2006; Müller et al., 2012b) and does not include the overburden of the material exterior to the shock. A comparison between the two is given in Fig. 11. There, we compute the diagnostic energy as the integral of the total energy density $e$ minus the zero-point energy $e_{0}$, $e_{\rm tot}=e-e_{0}$, either integrated over regions where $e_{\rm tot}\geq 0$, or over $100\ {\rm km}\leq r\leq 20,\!000\ {\rm km}$. For clarity, we did not include the binding energy of the envelope when computing $E_{\rm tot}$ in this plot, since its inclusion in the diagnostic explosion energy would be inconsistent. Beside this difference, the $r\geq 100\ {\rm km}$ diagnostic energies in Fig. 11 are identical to the estimated explosion energies in Fig. 10. Obviously, the diagnostic energy integrated only over $e_{\rm tot}\geq 0$ or $r\geq 100\ {\rm km}$ should converge to the same value after a sufficiently long time. This is indeed the case for most of our models, especially the ECSN/ECSN-like explosions, where the explosion is close to spherically symmetric. However, significant differences persist until the end of our simulations for some progenitors. For example, in the 11.0-ManyBody run, the shock starts expanding ${\sim}200\ {\rm ms}$ after bounce, and some material becomes unbound. However, the energy behind the shock becomes sufficient to overcome the overburden only at later times, when the PNS wind becomes violent enough to create high-entropy bubbles behind the shock and the initially bound material at $r\geq 100\ {\rm km}$ has accreted or has become unbound.

5. Neutrino Radiation

We extract the properties of the neutrino radiation on a sphere placed at $10,\!000\ {\rm km}$ from the center. Angle-averaged neutrino luminosities and rms neutrino energies are shown in Fig. 12. These are shown in the lab frame at infinity. Note that Fornax evolves the fluid-frame neutrino-radiation moments and does not output the Eddington factor used for the evolution. For this reason, we convert the code output to the lab frame under the simplifying assumption of a forward-peaked radial neutrino distribution function. This assumption is not valid at the time the shock crosses 10,000 km and results in small jumps which are particularly evident in the rms energies of heavy-lepton neutrinos for the n8.8 progenitor.

We find the ManyBody setup to result in slightly higher neutrino luminosities and average energies. The heavy-lepton neutrino luminosities are the most clearly affected and increase by ${\sim}10\%$, since their opacity is dominated by neutrino-nucleon scattering. However, the average energies for all neutrino species increase because of the accelerated contraction rate of the PNS with the ManyBody setup (see Sec. 6).

We quantify the degree of coupling between the neutrino radiation and the accretion flow in terms of the heating efficiency parameter $\eta$, defined as the ratio between the heating rate by neutrinos in the gain region, i.e., the region bounded by the PNS and the shock with positive net neutrino heating, and the sum of the $\nu_{e}$ and $\bar{\nu}_{e}$ luminosities at infinity (e.g., Marek & Janka, 2009; Müller et al., 2012b, a). We show this quantity as a function of time in Fig. 13 for 1D and 2D models with the Baseline and the ManyBody setups. We recall that the same analysis was performed by Müller et al. (2012a) for the u8.1 progenitor. We find good agreement with their heating efficiency and neutrino luminosities. For the other models, the overall trend in $\eta$ is that progenitors with larger accretion rates also show larger heating efficiencies.

Spherically-symmetric (1D) simulations have significantly smaller heating efficiencies. The increased heating efficiency in 2D is in part due to the longer dwelling time of material in the gain region (Burrows et al., 1995; Murphy & Burrows, 2008; Dolence et al., 2013) or, equivalently (Müller et al., 2012b) to the growth of the mass in the gain region.

The many-body corrections implemented in the ManyBody setup result in a slight increase of the heating efficiency. At least at early times, before the evolutionary paths of the Baseline and ManyBody simulations start to diverge, this improvement can be attributed to the hardening of the neutrino spectra with the ManyBody setup, which results in a more tight coupling with the material. After the explosions set in, the differences between the Baseline and ManyBody efficiencies are in good part due to the fact that the ManyBody explosions are more spherical and entrain more mass.

6. Protoneutron Stars

As is commonly done in the CCSN-mechanism literature, we define as PNS radius the radius at which the angle-averaged density is $10^{11}\ {\rm g}\ {\rm cm}^{-3}$. We monitor PNS radii and the baryonic mass they enclose to estimate the final remnant radius. These quantities are shown in Figs. 14 and 15. PNS masses and accretion rates at the end of our simulations are also given in Tab. 3. There we also quote the corresponding gravitational mass for a cold, deleptonized neutron star (NS) estimated using the approximate fit of Timmes et al. (1996).

The PNS masses, as was the case for the explosion energy, are still not converged for the 11.0-$M_{\odot}$ progenitor at the end of our simulations. The PNS mass for the 10.0-$M_{\odot}$ model is obviously converged only with the ManyBody setup with perturbations, which explodes, while black-hole formation would appear inevitable for the other setups.

Notwithstanding these caveats, we find that all our progenitors produce PNSs with gravitational masses below $1.4\ M_{\odot}$. The n8.8 progenitor produces PNSs with gravitational masses as low as $1.188\ M_{\odot}$, suggesting that an ECSN/ECSN-like explosion might be able to explain the origin of the low-mass companion in the double NS system J0453+1559. This has been recently measured, using the advance of periastron and the Shapiro delay, to have a mass of $1.174\pm 0.004\ M_{\odot}$ (Martinez et al., 2015). This scenario is also plausible in light of the study by Tauris et al. (2015), who showed that an ECSN is a possible outcome of the evolution of ultra-stripped metal cores in tight binaries, such as those producing relativistic double-NS systems like J0453+1559.

We caution the reader that previous studies reported somewhat larger PNS masses for the n8.8 progenitor. The Garching group reported a final (baryonic) PNS mass of $1.366\ M_{\odot}$ (Hüdepohl et al., 2010), while Fischer et al. (2010) reported a final PNS mass of $1.347\ M_{\odot}$. The origin for our smaller PNS masses is probably due to our neglecting of electron capture on heavy nuclei and nuclear burning during infall for this model. Both could slightly increase the PNS mass. As argued by Burrows & Lattimer (1985), nuclear burning in the supersonically infalling material will accelerate the collapse. Electron capture on heavy nuclei will decrease the pressure support in the core and further accelerate the collapse.

We find the PNS radii (Fig. 15) to follow tracks that are largely independent of the progenitor or the PNS mass, as do Bruenn et al. (2016) and Summa et al. (2016). The reason is that the density drops sharply at the surface of the PNS so that the ambient pressure has a negligible influence on the structure of the central object. This is determined by the competition between its internal pressure and gravity. Instead, the radii are sensitive to changes in the microphysical treatment, which determines the rate at which the PNS deleptonizes and loses thermal support, and to the dimensionality (1D vs. 2D). The impact of the microphysics is easily understood from the fact that the contraction of the PNS is mostly set by the rate of deleptonization and core cooling. These in turn depend in the first second after bounce on the neutrino opacity of matter at densities between $10^{11}$ and $10^{13}\ {\rm g}\ {\rm cm}^{-3}$. For instance, the many-body corrections included in the ManyBody setup result in a faster deleptonization and contraction of the PNS.

The reason for the faster PNS contraction in 1D is more easily understood considering the n8.8 progenitor with the Baseline setup, which explodes both in 1D and in 2D. Its PNS evolution is shown in Fig. 16. In the first few hundreds of milliseconds after bounce, the convection inside the PNS is buried deep below the surface and its impact on supernova evolution is limited, as has been documented in detail by Buras et al. (2006) and Dessart et al. (2006). However, over timescales longer than those considered in either of those works, starting from ${\sim}0.2\ {\rm s}$ after bounce, the surface of the PNS shrinks to the point of entering in contact with the inner PNS convection, which then becomes dynamically important. This can be seen in Fig. 16, where the region affected by the inner convection is identifiable by its small radial entropy gradient with entropy per baryon evolving from ${\sim}4.5\ {\rm k_{B}}$ to ${\sim}3\ {\rm k_{B}}$ as the PNS cools down.

Starting from this moment, the 1D and 2D evolutions begin to diverge. In 1D, the neutrino cooling of the surface is not compensated by convection and leads to an increasingly steep entropy inversion. The pressure support in the exterior layers of the PNS drops rapidly and leads to an increased compactness, with respect to the 2D evolution, of the regions with densities between $10^{11}$ and $10^{13}\ {\rm g}\ {\rm cm}^{-3}$. These regions are instead inflated in 2D by the deposition of entropy and lepton number due to convective transport. The structure of the layers below the inner convection region is also affected, with the core of the PNS reaching higher densities and compactness in 1D.

PNS convection also leaves a strong imprint on the neutrino luminosity, which is boosted by up to a factor ${\sim}2$ at late times ($\gtrsim 0.5\ {\rm s}$), as can be seen from Fig. 17. This seems to be the main reason for the enhanced growth of the explosion energy for the n8.8, u8.1, and z9.6 progenitors in 2D at late times. While this accounts only for a relatively small fraction of the explosion energy (see Fig. 10), the role of PNS convection might be more important for more massive progenitors that explode later in multi-dimensional simulations.

We remark that O’Connor & Couch (2015) also reported modest increases in the heavy-lepton neutrino luminosities due to PNS convection. However, since they considered models that did not explode in 1D, they might have underestimated the effect of convection, since the $\bar{\nu}_{e}$ and, in particular, the $\nu_{e}$ luminosities are significantly affected by accretion. Indeed, while the z9.6 and u8.1 progenitor evolutions (the former only with the ManyBody setup) show very similar differences between 1D and 2D as does the n8.8, this is not the case for the 9.0-, 10.0-, and 11.0-$M_{\odot}$ progenitors, which show more similar luminosities in 1D and 2D.

The n8.8 model was considered in 1D by Fischer et al. (2010), who found essentially the same luminosity as in our 1D n8.8-Baseline model ($L_{\nu_{e}}\simeq\ 6\cdot 10^{52}\ {\rm erg}\ {\rm s}^{-1}$ at $0.6\ {\rm s}$ after bounce). This is, however, a factor ${\sim}2$ smaller than in our 2D calculations for the n8.8. Müller & Janka (2014) also considered the z9.6 progenitors over a long timescale and found luminosities very close to ours with the Baseline setup. Their luminosity was $L_{\nu_{e}}\sim 10^{52}\ {\rm erg}\ {\rm s}^{-1}$ at $0.6\ {\rm s}$ after bounce, the same value we also find (Fig. 12), but they did not present a comparison with the corresponding 1D evolution. These considerations are all additional indirect confirmation that the impact of PNS convection has been underestimated.

7. Conclusions

We have revisited the explosion of low-mass iron-core SNe and O-Ne-Mg gravitational-collapse SNe with a new set of neutrino-radiation hydrodynamics simulations in 1D (spherical symmetry) and in 2D (axial symmetry). Our simulations included the effects of general relativity in an approximate way and state-of-the-art multidimensional neutrino transport and weak reactions. Of the six progenitors we have considered, one, the n8.8 from Nomoto (1984, 1987), is the prototype of an ECSN. Two, the $10^{-4}\ Z_{\odot}$ and zero metallicity u8.1 and z9.6 from Heger (private communication), have iron cores, but a structure similar to that of the n8.8. The 9.0-, 10.0-, and 11.0-$M_{\odot}$ solar-metallicity progenitors from Sukhbold et al. (2016) share some similarities with the n8.8, but are overall closer to the “canonical” CCSN progenitors considered in in the CCSN-mechanism literature.

As in previous studies (Kitaura et al., 2006; Janka et al., 2008; Burrows et al., 2007; Fischer et al., 2010; Müller et al., 2012a; Janka et al., 2012; Müller et al., 2013; Müller & Janka, 2014; Melson et al., 2015b; Wanajo et al., 2017), we find that the n8.8 and z9.6 progenitors typically explode easily, even in 1D. The u8.1 progenitor is close to the threshold for explosion in 1D and is successful when many-body corrections to neutral current reactions are included, as in our ManyBody setup. On the other hand, the low-mass, but solar-metallicity, progenitors with iron cores from Sukhbold et al. (2016) do not explode in 1D and, in some cases, like the 10.0-$M_{\odot}$, also fail to explode in 2D. Our results show that solar-metallicity iron-core SNe do not explode in 1D and are not even necessarily easier to explode than higher-mass stars.

The failure to explode of the 10.0-$M_{\odot}$ progenitor from Sukhbold et al. (2016) reported here is surprising in the light of the successful explosion of other progenitors, like the 9.0- and 11.0-$M_{\odot}$ progenitors from Sukhbold et al. (2016). Our findings are in tension with explodability criteria related to the ZAMS mass (Heger et al., 2003) or to the progenitor compactness and related parameters (O’Connor & Ott, 2011, 2013; Ugliano et al., 2016; Nakamura et al., 2015; Pejcha & Thompson, 2015; Ertl et al., 2016). Other circumstantial evidence of limitations in the existing explodability criteria is the order in which explosions develop in Summa et al. (2016) and the results of O’Connor & Couch (2015). The latter considered the 12-, 15-, 20-, and 25-$M_{\odot}$ progenitors from Woosley & Heger (2007) and found explosions in approximate GR for all progenitors except the $12\ M_{\odot}$. Taken together, all of these results suggest that, while some properties of the explosions are correlated with the compactness of the progenitor (e.g., O’Connor & Ott, 2013; Nakamura et al., 2015), the explodability is not. Whether an explosion is successful or not depends on a competition between accretion and neutrino heating (Burrows & Goshy, 1993; Janka, 2000; Suwa et al., 2016; Murphy & Dolence, 2017; Gabay et al., 2015) which, in our opinion, has yet to be expressed in terms of the progenitor properties in a satisfactory way.

We have systematically studied the effect of perturbations and of changes in the treatment of neutrino-matter interactions, with emphasis on the impact of the many-body corrections to the neutrino-nucleon scattering cross section derived by Horowitz et al. (2017). We have found that relatively small changes are amplified by the proximity to the threshold for explosion and can lead to qualitatively different outcomes. For instance, the 10.0-$M_{\odot}$ model turns from a dud into an explosion with the inclusion of perturbations in combination with many-body corrections. This sensitivity to initial conditions and/or physics setup applies also in 1D to those progenitors that are sufficiently close to the threshold for explosion. For example, the u8.1 fails to explode in 1D without the inclusion of many-body corrections, but succeeds when those are included. The reason for the diverging outcomes with different microphysical suites is easily understood in terms of the neutrino-radiation intensity and hardness, which directly translate into the efficiency of the energy deposition by neutrinos.

We have estimated explosion energies by following the development of the explosions over long timescales and until the shock has reached 19,000 km in 2D simulations. While the explosion energy for the 11.0-$M_{\odot}$ progenitor is still far from saturated, for the others we have found saturated explosion energies of the order of a tenth of a Bethe. These values are in the expected range for low-mass progenitors (Utrobin & Chugai, 2013; Spiro et al., 2014; Sukhbold et al., 2016). We remark that, while the ECSN/ECSN-like explosions are nearly-spherical and we do not expect their explosion energies to change significantly in 3D (Melson et al., 2015b), it is likely that the asymmetric explosions we observe for the 9.0-, 10.0-, and 11.0-$M_{\odot}$ progenitors will be quantitatively different in 3D (Müller, 2015).

For the progenitors that have been evolved in the past by other groups (the n8.8, u8.1, and z9.6) we obtain explosion energies in 2D that are typically ${\sim}50\%$ larger than those of the Garching group (Janka et al., 2008; Melson et al., 2015b; Wanajo et al., 2017). On the other hand, we find good agreement with the Garching results for the explosion energy of the z9.6 progenitor in 1D, which suggests that these difference might be ascribed to multi-dimensional effects, such as the handling of convection inside the PNS and behind the shock and the problematic ray-by-ray approximation for lateral neutrino transport (Skinner et al., 2016). We also remark that our calculations did not include nuclear burning and instead assumed nuclear statistical equilibrium, which could also explain some of the differences. We did not attempt to systematically investigate the reasons for the discrepancies, but suggest that the use of models that explode in self-consistent 1D simulations, in the context of studies on the impact of the microphysics on the explosion mechanism, or within a renewed effort to cross-validate CCSN codes, appears promising.

For exploding models, we have found final PNS masses that have reached saturation within the simulation time, with the exception of the 11.0-$M_{\odot}$ progenitor. We find that ECSNe can explain the low-mass tail of the observed NS mass distribution. The n8.8 progenitor with the ManyBody setup leaves behind a NS with baryonic (gravitational) mass of $1.294\ M_{\odot}$ ($1.188\ M_{\odot}$), very close to the lowest accurately measured NS gravitational mass of $1.174\pm 0.004\ M_{\odot}$ (Martinez et al., 2015).

We studied the evolution of the PNS, focusing on those progenitors that explode both in 1D and 2D. These have nearly identical boundary conditions, allowing us to quantify the role of multi-dimensional effects on the long-term evolution of the PNS. We have found that the PNS contraction rate slows down significantly in 2D compared to 1D, starting from ${\sim}0.2\ {\rm s}$ after bounce. At this time, the PNS surface has contracted sufficiently to enter in contact with the inner PNS convection region. The transport of lepton number and thermal energy by the PNS convection then inflates the region with densities between $10^{11}$ and $10^{13}\ {\rm g}\ {\rm cm}^{-3}$, causing a decrease of the contraction rate.

We have also found PNS convection to be responsible for a boost of the neutrino luminosities for all species by up to a factor ${\sim}2$ at late times $\gtrsim 0.5\ {\rm s}$. This contributed only a $\lesssim 10\%$ increase to the explosion energy for the ECSN/ECSN-like progenitor, for which this effect is more easily quantifiable. However, PNS convection is likely to be more important for massive progenitors that explode late, when the PNS surface has already receded sufficiently close to the PNS convection region. Our results, together with pieces of evidence from Fischer et al. (2010) and Müller & Janka (2014), strongly suggest that the impact of PNS convection has been underestimated in the past (Buras et al., 2006; Dessart et al., 2006). Our findings provide an additional reason, beside the need to account for continued accretion at late times (Müller & Janka, 2014), for the importance of multi-D simulations in the modeling of the early neutrino signal from cooling PNSs (e.g., Fischer et al., 2010; Hüdepohl et al., 2010; Roberts et al., 2012b; Roberts, 2012; Nakazato et al., 2013; Roberts & Reddy, 2016). At the same time, we caution the reader that our simulations did not include in-medium modifications of charged-current reactions, which will also affect the quantitative properties of the PNS neutrino-cooling light curve, especially after the first second (Burrows & Sawyer, 1999; Martinez-Pinedo et al., 2012; Horowitz et al., 2012; Roberts et al., 2012a).

The main limitation of this work is the assumption of axisymmetry. This has been necessary given the large computational cost of 3D simulations with state-of-the-art microphysics, which prevents a systematic study with multiple progenitors, as is the present one. This is attested by the scarcity of 3D simulations with full microphysics (Hanke et al., 2013; Tamborra et al., 2014; Melson et al., 2015b, a; Lentz et al., 2015; Summa et al., 2017). However, moving to 3D will ultimately be required and will be a goal of our future work.

Acknowledgments

The authors acknowledge Chuck Horowitz, Evan O’Connor, and Todd Thompson for productive conversations concerning, insight into, and help with the microphysics, and Ernazar Abdikamalov, Sean M. Couch, Luke F. Roberts and Christian D. Ott for fruitful discussion on the nature of core-collapse supernovae. Support was provided by the Max-Planck/Princeton Center (MPPC) for Plasma Physics (NSF PHY-1523261). DR gratefully acknowledges support from the Schmidt Fellowship. AB acknowledges support from the NSF under award number AST-1714267. JD acknowledges support from a Laboratory Directed Research and Development Early Career Research award at LANL. The authors employed computational resources provided by the TIGRESS high performance computer center at Princeton University, which is jointly supported by the Princeton Institute for Computational Science and Engineering (PICSciE) and the Princeton University Office of Information Technology and by the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the US Department of Energy (DOE) under contract DE-AC03-76SF00098. The authors express their gratitude to Ted Barnes of the DOE Office of Nuclear Physics for facilitating their use of NERSC. This paper has been assigned a LANL preprint # LA-UR-17-20973.