The Progenitor Dependence of Core-collapse Supernovae from Three-dimensional Simulations with Progenitor Models of 12–40 M_⊙

Core-collapse supernovae (CCSNe) are the birth places of neutron stars and black holes. They liberate the ashes of stellar evolution, seeding the interstellar gas with the elements from which planets form and life is made. They feed back on star formation and regulate galaxy gas budgets. Yet, despite their importance for much of astrophysics, our understanding of the CCSN explosion mechanism, and its dependence on progenitor star properties, is woefully incomplete.

The CCSN problem (see, e.g., Janka et al. 2007 for an in-depth review) boils down to the total pressure behind the shock having to offset the accretion ram pressure of the outer core impinging on the shock. The hot accreting protoneutron star (PNS) formed at core bounce emits a huge flux of neutrinos of all species. The neutrino mechanism (Bethe & Wilson 1985) relies on a fraction of these neutrinos being reabsorbed in a "gain layer" below the shock. There, they heat the gas, increasing the thermal pressure. Stars are spherical from a distance and much of the early CCSN simulation work was conducted in spherical symmetry (1D). But 1D simulations fail to show explosions powered by the neutrino mechanism for all but the lowest-mass progenitors (M ≲ 10 M_⊙; e.g., Kitaura et al. 2006; Radice et al. 2017). Neutrino heating is strongest at the base of the gain layer and it establishes a radially decreasing gradient in entropy. The first axisymmetric (2D) simulations showed that turbulent convection driven by this gradient plays a crucial role in reviving shock expansion (Herant et al. 1994; Burrows et al. 1995; see Couch & Ott 2015 for the role of turbulence). 2D simulations also showed that a non-spherical instability of the standing-accretion shock (SASI; Blondin et al. 2003; Foglizzo et al. 2006) can also help revive the shock.

The recent availability of petascale supercomputers has enabled the first detailed 3D CCSN simulations (Hanke et al. 2013; Takiwaki et al. 2014, 2016; Tamborra et al. 2014; Lentz et al. 2015; Melson et al. 2015a, 2015b, Roberts et al. 2016; Müller et al. 2017; Chan et al. 2018; Summa et al. 2018). Comparisons with 2D simulations have shown that 3D is essential for understanding CCSNe, their explosion mechanism, and for predicting their multi-messenger observables (see Couch & Ott 2015; Janka et al. 2016).

In this Letter, we present a first study of the progenitor-star dependence of neutrino-driven CCSNe in 3D, covering zero-age main-sequence masses from 12 M_⊙ to 40 M_⊙. All progenitors, except for the 12 M_⊙ star, see shock runaway by 500 ms after bounce, but in remarkably distinct ways, depending sensitively on their precollapse structure.

We draw 1D progenitors of 12, 15, 20, and 40 M_⊙ from the set of Woosley & Heger (2007, hereafter WH07). In addition, we use the 27 M_⊙ model of Woosley et al. (2002, hereafter WHW02) that was simulated in 3D by Hanke et al. (2013) and Roberts et al. (2016). We plot the progenitor density profiles in Figure 1.

Zoom In Zoom Out Reset image size

We simulate core collapse in 1D using GR1D (O'Connor & Ott 2013; O'Connor 2015) and map to 3D at 20 ms after bounce for all WH07 progenitors. The 27 M_⊙ model is mapped at 38 ms after bounce due to a transposed-digits error of the lead author. We carry out the 3D simulations with the open-source 3D general-relativistic (GR) multi-group radiation-hydrodynamics CCSN code Zelmani (Roberts et al. 2016). It is based on the Einstein Toolkit (Löffler et al. 2012; Mösta et al. 2014). Neutrino transport is handled in the GR M1 multi-group approximation (Shibata et al. 2011). We use three neutrinos species (ν_e, ${\bar{\nu }}_{e}$, and ${\nu }_{x}=[{\nu }_{\mu },{\bar{\nu }}_{\mu },{\nu }_{\tau },{\bar{\nu }}_{\tau }]$) and 12 energy groups, spaced logarithmically with bin centers between 1 and 248 MeV. In the 3D simulations, we employ the subset of Bruenn (1985) neutrino opacities used in O'Connor & Ott (2013), but leave out velocity dependence and inelastic scattering processes. Velocity dependence and inelastic neutrino–electron scattering are included in the 1D collapse simulations as described in O'Connor (2015). All simulations employ the SFHo equation of state, which is tuned to fit astrophysical and experimental constraints (Steiner et al. 2013). Upon mapping from 1D, we observe a short transient in all neutrino quantities for about a light-crossing time until the radiation field has reached quasi-steady state and settles within ∼5% of the 1D results. Similarly, there is a small re-adjustment in the hydrodynamics as the 1D structure relaxes to the 3D Cartesian grid.

The 3D simulations use eight levels of Cartesian adaptive mesh refinement, resolving the PNS with 370 m resolution and the postshock region with 1.5 km before shock expansion (see Abdikamalov et al. 2015; Roberts et al. 2016 for resolution studies with our code). After the shock has expanded to radii ≳300 km, we regrid to 3 km resolution for the shocked region. We find that once the shock has begun its runaway, changes in the postshock resolution have only small effects. Table 1 summarizes key model properties. All times in this Letter are measured relative to core bounce of each model.

Table 1.  Model Summary

Model	ξ1.75	Mic,b (M⊙)	tf − tb (ms)	$\langle {R}_{\mathrm{shock},{\rm{f}}}\rangle \,(\mathrm{km})$
s12WH07	0.235	0.583	527	123
s15WH07	0.580	0.576	597	526
s20WH07	0.944	0.577	384	523
s27WHW02	0.783	0.573	392	482
s40WH07	1.328	0.562	323	614

Note. ξ_1.75 is the core compactness (O'Connor & Ott 2011) measured at bounce at a mass coordinate of 1.75 M_⊙. M_ic,b is the mass of the homologous core at bounce. t_f − t_b is the final simulation time relative to core bounce and $\langle {R}_{\mathrm{shock},{\rm{f}}}\rangle$ is the final average shock radius.

Download table as:  ASCIITypeset image

Core bounce occurs when the inner core reaches nuclear density and the repulsive short-range component of the nuclear force stabilizes its collapse. With little variation between progenitors, the CCSN shock is launched from a mass coordinate of 0.56–0.58 M_⊙ (Table 1; this is expected, see, e.g., Janka et al. 2012). The shock first expands rapidly, but quickly weakens due to the dissociation of heavy nuclei in accreting outer core material and neutrino losses. It succumbs to the accretion ram pressure and stalls at a radius of ∼150 km. At this point, differences in progenitor structure begin to matter.

In the bottom left panel of Figure 2, we plot the time evolution of the mass accretion rate $\dot{M}$ in all progenitors. Within tens of milliseconds of bounce, the entire iron core has accreted in all models. Comparing with Figure 1, we see that the subsequent $\dot{M}$ is determined by the progenitor density profile in the overlying Si and Si–O shells at mass coordinate M ≳ 1.3–1.9 M_⊙.

Zoom In Zoom Out Reset image size

One expects the accretion rate to have multiple important roles, some of which may counteract one another. First, it sets the accretion ram pressure ${P}_{\mathrm{ram}}=\rho {v}^{2}\propto \dot{M}{M}_{\mathrm{PNS}}^{1/2}{r}_{s}^{-5/2}$ (the second relationship results from assuming the accreted material is in free-fall from a large radius), which keeps the stalled shock from expanding. Second, it regulates the accretion luminosity ${L}_{\mathrm{acc}}\propto \dot{M}{M}_{\mathrm{PNS}}{R}_{\mathrm{PNS}}^{-1}$, which is the dominant source of ν_e and ${\bar{\nu }}_{e}$ luminosity providing energy to the shock. The hierarchy of ν_e and ${\bar{\nu }}_{e}$ luminosities between the different progenitor models shown in the center panels of Figure 2 directly reflects the $\dot{M}$ order. Finally, the integrated accretion rate determines mass and radius evolution of the PNS. To within an order of magnitude, the mean neutrino energies are proportional to the temperature at the surface of the PNS. A virial argument suggests $\langle {\epsilon }_{{\nu }_{e}}\rangle \propto {({M}_{\mathrm{PNS}}{R}_{\mathrm{PNS}}^{-1})}^{\alpha }$, which will depend on the integrated accretion rate.⁸ At early times, ${M}_{\mathrm{PNS}}{R}_{\mathrm{PNS}}^{-1}$ is similar among all progenitors due to the universal structure of the collapsed cores (Figure 3) and therefore the mean neutrino energies are similar (Figure 2), but they start to deviate with time due to differing accretion histories. Because of all that, one expects the accretion history to be a determining factor in the CCSN evolution of a given progenitor.

Zoom In Zoom Out Reset image size

The impact of progenitor structure on the shock evolution toward explosion is depicted by the top left panel of Figure 2. We see three qualitatively different evolution modes:

• (1)  
  The 40 M_⊙ progenitor has the highest postbounce $\dot{M}$, translating to the highest neutrino luminosities. Its density profile (Figure 1) is shallow and smooth and there are no quick drops in $\dot{M}$. This model's shock begins to deviate substantially from spherical symmetry at ∼100 ms after bounce and shock runaway ensues at around 200 ms.
• (2)  
  The 20 M_⊙ and 27 M_⊙ models have lower postbounce $\dot{M}$, but their density profiles have a steep discontinuity at the Si/Si–O shell interface⁹ (see Figure 1). In both models, it is the drop in ram pressure due to the rapidly decreasing $\dot{M}$ that triggers shock runaway ∼170–200 ms after bounce.
• (3)  
  In the 12 M_⊙ and 15 M_⊙ models with their moderate $\dot{M}$ and low L_ν, the shock recedes to radii around 100 km. The accretion rate gradually decreases, and so do the ν_e and ${\bar{\nu }}_{e}$ luminosities (center panels of Figure 2), while the mean neutrino energies increase due to the increasing compactness of the PNS (bottom right panel of Figure 3). Both models experiences SASI. Eventually, more than 500 ms after bounce, shock runaway occurs in the 15 M_⊙ model. The 12 M_⊙ model does not experience shock runaway by the end of our simulation, but it still has the potential to resume expansion at a later time.

O'Connor & Couch (2015) and Summa et al. (2016) found similar evolutions to modes (2) and (3) in 2D simulations.

In Figure 3, we present diagnostics that help understand the three evolution modes. Shock expansion is facilitated by increases in thermal and turbulent pressure that offset the accretion ram pressure (e.g., Couch & Ott 2015). Stronger neutrino heating means more thermal pressure and stronger driving of turbulence. The neutrino heating rate scales approximately as ${\dot{Q}}_{\mathrm{heat}}\propto (\langle {\epsilon }_{{\nu }_{e}}^{2}\rangle {L}_{{\nu }_{e}}+\langle {\epsilon }_{{\bar{\nu }}_{e}}^{2}\rangle {L}_{{\bar{\nu }}_{e}}){R}_{\mathrm{gain}}^{-2}{M}_{\mathrm{gain}},$ where M_gain and R_gain are the mass contained in the gain region and the gain radius, respectively (Janka 2001; Summa et al. 2016). Therefore, the hierarchy of heating rates among the models mirrors their luminosity hierarchy (Figure 3). Assuming that the majority of the ν_e and ${\bar{\nu }}_{e}$ luminosity is powered by accretion, one finds ${\dot{Q}}_{\mathrm{heat}}\propto \dot{M}{({M}_{\mathrm{PNS}}{R}_{\mathrm{PNS}}^{-1})}^{1+2\alpha }{R}_{\mathrm{gain}}^{-2}{M}_{\mathrm{gain}}$, which implies greater heating for a higher accretion rate and a more compact PNS for a fixed gain region size and mass. Interestingly, at early times (≲80–100 ms), the heating efficiency $\eta ={\dot{Q}}_{\mathrm{net}}{({L}_{{\nu }_{e}}+{L}_{{\bar{\nu }}_{e}})}^{-1}$ (where ${\dot{Q}}_{\mathrm{net}}$ is the net heating rate; heating minus cooling) is independent of progenitor. Since the mean neutrino energies are very similar at early times, this implies that ${M}_{\mathrm{gain}}{R}_{\mathrm{gain}}^{-2}$ is similar for all of the models even though they have different masses in the gain region (see the top center panel of Figure 3).

Neutrino-driven convection begins to grow at ∼80–100 ms in all models, as can be seen from the top right panel of Figure 3, showing the radial and nonradial specific turbulent kinetic energy in the gain region. It is fully developed at 200 ms (Figure 4).

Zoom In Zoom Out Reset image size

The specific turbulent energy in the gain region is very similar in all models and grows with time, although convection starts somewhat later for s40WH07. When turbulent convection begins in each model, it briefly reverses the decline of M_gain. In most models, the decline in M_gain continues, but in s40WH07, the onset of neutrino-driven turbulent convection stabilizes M_gain, as shown by the top center panel of Figure 3. At constant M_gain and near-constant luminosity, the increase in neutrino heating in s40WH07 (top left panel of Figure 3) follows from the steep postbounce increase of the neutrino energies shown in Figure 2.

The onset of neutrino-driven convection in model s40WH07 around 100 ms after bounce is also reflected in the departure of its shock from spherical symmetry and the expansion of its maximum shock radius (Figure 2). Conditions for global shock runaway become gradually more favorable. This can be seen by comparing the timescale ${\tau }_{\mathrm{adv}}\approx {M}_{\mathrm{gain}}{\dot{M}}^{-1}$ it takes for material to advect through the gain layer with ${\tau }_{\mathrm{heat}}\approx | {E}_{\mathrm{gain}}| {\dot{Q}}_{\mathrm{net}}^{-1}$, the timescale for neutrino heating (Janka 2001; Thompson et al. 2005; we follow Müller et al. 2012 for implementation details). If τ_adv/τ_heat ≳ 1, it is said that conditions favor shock runaway. From the bottom center panel of Figure 3, we see that s40WH07 crosses τ_adv/τ_heat = 1 at ∼200 ms after bounce. This is at about the time its shock runs away globally. However, τ_adv/τ_heat ≳ 1 seems to be more of a diagnostic rather than a critical condition for explosion. Instead, what appears to set model s40WH07 on a positive track toward shock runaway is the stabilization of the mass in its gain layer. This occurs when neutrino-driven turbulence sets in, about 100 ms before s40WH07 reaches τ_adv/τ_heat ∼ 1.

The diagnostics for models s12WH07, s15WH07, s20WH07, and s27WHW02 shown in Figure 3 are qualitatively very similar until ∼150–180 ms after bounce. Neutrino heating scales roughly with the luminosities, which in turn scale with the accretion rates (Figure 2). M_gain decreases with the accretion rate and the heating efficiencies show only moderate dependence on progenitor model.

Models that fall into mode (2) defined above, s20WH07 and s27WHW02, depart from this trend when their Si/Si–O interface reaches the shock. In s27WHW02, this occurs at ∼180 ms, leading to the steep drop in $\dot{M}$ seen in the bottom left panel of Figure 2. Correspondingly, the ram pressure lid on the shock is lifted, M_gain stabilizes, τ_adv/τ_heat (Figure 3, bottom center panel) jumps above 1, and global shock expansion ensues. An important aspect is that for one advection time of τ_adv ∼ 10–15 ms, the neutrino luminosity and heating rate remain approximately constant. They drop only once the lower $\dot{M}$, due now to both the density drop and the shock expansion, has propagated to the PNS edge. This helps power shock expansion. A drop of $\dot{M}$ at near-constant luminosity was identified, e.g., by Suwa et al. (2016), as conducive to explosion. Model s20WH07 mirrors what we find for s27WHW02.

Models s12WH07 and s15WH07 also mirror each other for most of their evolution (see Figure 3). Their accretion rates are moderate and neutrino heating and, consequently, neutrino-driven turbulence is not strong enough to keep the shock on a positive trajectory. In the following, we focus on s15WH07.

After ∼150 ms, the shock begins to recede despite near-constant heating rates. τ_adv/τ_heat hovers around 0.5. The average shock radius settles around 100 km and remains there for hundreds of milliseconds, as shown by the top left panel of Figure 2. We note that the shock departs substantially from spherical symmetry with minimum and maximum shock radii differing by ∼40 km in s15WH07 at 200–500 ms after bounce.

Shock recession leads to short advection times and conditions favorable for the growth of SASI (e.g., Scheck et al. 2008; Fernández 2015). SASI with a substantial spiral mode (Fernández 2010) develops in s15WH07 and s12WH07. In Figure 4 (second row, first column), we show a snapshot of the lopsided specific entropy distribution in the x–y plane of model s15WH07 at 393 ms after bounce, which shows a SASI-deformed shock. In the right panel of Figure 5, we plot the normalized ℓ = 1, m = {−1, 0, 1} mode amplitudes of a spherical harmonics expansion of the shock front. The mode amplitudes grow 100–200 ms after bounce when neutrino-driven convection is present (see top and second row, first column of Figure 4). They become clearly oscillatory and persistent once SASI develops, ∼250 ms after bounce.

Zoom In Zoom Out Reset image size

With the help of SASI and a growing specific turbulent energy in the gain region, neutrino heating gradually brings s15WH07 back onto a positive track toward shock revival. Figure 3 shows that the mass M_gain contained in the gain region stabilizes around 300 ms after bounce. While the neutrino luminosity continues its slow decline (Figure 2), heating rate and heating efficiency both increase since M_gain is stabilized and neutrino energies increase. The timescales ratio τ_adv/τ_heat follows and slowly reaches unity around 400 ms after bounce. Yet, shock expansion does not follow immediately and s15WH07 straddles the threshold of shock runaway for another ∼100 ms until heating and turbulence are strong enough to overcome the slowly decreasing ram pressure.

In contrast to s15WH07, model s12WH07 does not experience shock runaway within the simulation time. However, at the end of the simulation, s12WH07 has τ_adv/τ_heat > 1 and the mass in the gain region has stabilized, both of which suggest it is on the path to shock revival.

Once the average shock radius reaches ∼300 km, shock expansion rapidly accelerates in all models. In the third row of Figure 4, we show x–y slices of the specific entropy at the time the average shock radius reaches 300 km in models s20WH07, s27WHW02, and s40WH07. The expanding shock already exhibits large asymmetry in these models, which only grows as the expansion accelerates. This can be appreciated from the third and forth rows of Figure 4.

The asymmetry of the expanding shock can be understood more quantitatively by considering Figure 5. There, we plot the normalized rms mode amplitudes (summed over all m) for ℓ = 1 (top panel) and ℓ = 2 (bottom panel). As seen in previous exploding 3D simulations (e.g., Lentz et al. 2015; Melson et al. 2015a; Roberts et al. 2016; Müller et al. 2017), all models develop a strong, growing, and dominant ℓ = 1 asymmetry. There is also power in ℓ = 2 and in higher modes (not shown), reflecting the complex morphology of the expanding shock shown in the two bottom rows of Figure 4.

In the top panel of Figure 6, we plot the diagnostic energy (Equation (2) of Müller et al. 2012) of unbound material in our models. All our models are in the rapid shock expansion phase at the end of their simulations and the diagnostic energy grows quickly. It does not factor in the "overburden" (Bruenn et al. 2016), the binding energy of the overlying stellar material. Using the PNS masses in Figure 6 as approximate mass cuts, we find, based on precollapse structures, overburdens of (17, 7.7, 6, 3.3, 1.3) × 10⁵⁰ erg for the (40, 27, 20, 15) M_⊙ models, respectively. Hence, the diagnostic energies still must grow before the overall energy budget becomes positive and typical CCSN energies are reached. The PNS, via its neutrino-driven wind, will continue to inject energy into the explosion for seconds. Luminosity from continuing accretion, facilitated by the highly asymmetric shock expansion, will also contribute. Bruenn et al. (2016; in 2D) and Müller et al. (2017; in 3D) have shown that in this way typical final explosion energies of 10⁵⁰–10⁵¹ erg can be reached by the neutrino mechanism.

Zoom In Zoom Out Reset image size

Model s40WH07, however, is a special case. Its PNS has a gravitational (baryonic) mass of ∼2.05 M_⊙ (∼2.13 M_⊙) at the end of the simulation and it is still accreting at a rate of ∼0.45 M_⊙ s⁻¹ due to the asymmetry of the explosion. The cold-NS maximum gravitational mass supported by the SFHo EOS is 2.06 M_⊙. The hot PNS in our non-exploding 1D simulation of this progenitor collapses to a black hole at a gravitational mass of 2.33 M_⊙. Hence, unless the accretion rate drops substantially, model s40WH07 will form a black hole within ∼600–700 ms from the end of our simulation. This will shut off energy injection into the expanding shock, resulting in an anemic explosion or complete failure (Chan et al. 2018).

The gravitational PNS mass in the other models has more or less leveled off at the end of the simulation due to the competition of neutrino cooling with moderate amounts of continuing accretion. Using E_bind ≈ 0.084 M_⊙ c² (M_grav/M_⊙)² (Lattimer & Prakash 2001) and the bound baryonic masses shown in Figure 6, we estimate (lower limit) final remnant NS masses of (1.59, 1.65, 1.54) M_⊙ for s15WH07, s20WH07, and s27WHW02, respectively.

Simple answers are not to be had in 3D CCSN theory. Our simulations suggest a complicated interplay of accretion rate, neutrino heating, and 3D fluid dynamics that determines the resulting CCSN dynamics and final outcome.

The three "CCSN evolution modes" we identify depend on progenitor structure as follows:

• (1)  
  Massive cores with high compactness (ξ_M = (M/M_⊙) (R[M]/1000 km)⁻¹, M = 1.75–2.5 M_⊙; O'Connor & Ott 2011) and without large density drops at shell interfaces develop early neutrino-driven, turbulence-facilitated shock runaway, but likely make black holes, with or without the shock exploding the star.
• (2)  
  Cores with intermediate compactness and with a substantial density drop at the Si/Si–O interface develop neutrino-driven, turbulence-facilitated explosions when this interface reaches the shock. They make relatively massive NSs with M ≳ 1.5 M_⊙.
• (3)  
  Cores with moderate to low compactness and without precipitous density drop have receding shocks that develop SASI and run away only at late times once SASI, neutrino heating, and turbulence have established favorable conditions. Due to the late explosions, the resulting NSs are also relatively massive (M ≳ 1.4–1.5 M_⊙).

Modes (2) and (3) were seen previously in the 2D simulations of O'Connor & Couch (2015) and Summa et al. (2016). Mode (1), for the most extreme progenitors like s40WH07, is new. Pan et al. (2017) recently simulated this progenitor in 2D with the same EOS, but did not find an explosion. Chan et al. (2018) simulated a different 40 M_⊙ progenitor in 3D and modified neutrino opacities to obtain an explosion. In our simulation of s40WH07, turbulence driven by neutrino heating is essential for creating conditions allowing shock runaway. The simulation of Pan et al. (2017) appears to have much weaker turbulence. The strength of CCSN turbulence is sensitive not only to neutrino heating, but also to the magnitude of perturbations that enter through the shock (e.g., Couch & Ott 2013; Müller & Janka 2015; Müller et al. 2017). Hence, the differences between our simulations and those of others could be due to the relatively large numerical perturbations imposed by our Cartesian grid (e.g., Ott et al. 2013). This could, perhaps, explain why we find explosions for s27WHW02 and s20WH07 that did not explode in the spherical-coordinate 3D simulations of Hanke et al. (2013) and Tamborra et al. (2014), and Melson et al. (2015a), respectively. Another piece of evidence for this argument is that our simulations of s20WH07 and s27WHW02 are at all times closer to explosion than their 2D counterparts in Summa et al. (2016), despite the fact that other studies have shown that 2D is more conducive to explosion than 3D due to its (unphysical) inverse turbulent cascade (e.g., Couch 2013; Couch & Ott 2015; Lentz et al. 2015). However, one should keep in mind that Hanke et al. (2013), Tamborra et al. (2014), Melson et al. (2015a), and Summa et al. (2016) used a different EOS, as well as different approximations to the neutrino transport and different neutrino opacities.

Presupernova stars in the wild have physical perturbations in their iron cores, Si, and Si–O shells. Determining their properties requires full 3D stellar evolution simulations of the final phase before core collapse (Couch et al. 2015; Müller et al. 2016; Cristini et al. 2017).

Our simulations show that the development of large-scale asymmetric explosions with dominant ℓ = 1 components is a generic outcome and independent of progenitor in the mass range considered here. While we do not investigate them here, NS and black hole birth kicks require such asymmetric mass ejection (e.g., Janka 2013; Müller et al. 2017).

Our three CCSN evolution modes produce black holes or massive NSs (M ≳ 1.4–1.5 M_⊙). Hence, there must be other CCSN modes for producing lower-mass NSs. Low-mass progenitors with O–Ne and the lowest-mass iron core progenitors (M ≲ 10 M_⊙) explode even in 1D due to their very steep density profiles (e.g., Kitaura et al. 2006; Melson et al. 2015b; Radice et al. 2017). They could be responsible for low-mass NSs. They have less mass to eject and their explosions are also likely to be less asymmetric, leading to low kick velocities.

There are many ingredients to the CCSN phenomenon. Here, we investigated for the first time the progenitor dependence in 3D. Others have recently investigated the role of rotation on the 3D neutrino mechanism (e.g., Takiwaki et al. 2016; Summa et al. 2018). More such studies are needed to also investigate magnetohydrodynamic effects, the impact of various neutrino transport approximations (see Richers et al. 2017 for recent progress there), differences in microphysics (EOS, neutrino interactions), and other numerical issues such as hydrodynamics methods, grid geometries, and resolution.

We acknowledge helpful discussions with D. Radice, H. Nagakura, Y. Suwa, K. Kiuchi, M. Shibata, J. Takeda, H.-T. Janka, R. Bollig, M. Obergaulinger, S. Couch, E. O'Connor, P. Mösta, and K. S. Thorne. This research was partially supported by NSF grants CAREER PHY-1151197, PHY-1404569, OAC-1550514, and AST-1333520, and the Sherman Fairchild Foundation. The simulations took over a year to complete and were carried out on ${ \mathcal O }$(10,000) compute cores of NSF/NCSA Blue Waters (PRAC ACI-1440083, OCI-0725070 and ACI-1238993, and Illinois allocation baov), of Edison at the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231, and of the Texas Advanced Computing Center Stampede and Stampede2 clusters under NSF XSEDE allocation TG-PHY100033. This article has been assigned Yukawa Institute for Theoretical Physics report No. YITP-17-122.