Hydrodynamical Neutron-star Kicks in Electron-capture Supernovae and Implications for the CRAB Supernova

Observations of young radio pulsars and measured orbital parameters of neutron stars (NSs) in close binary systems yield evidence that NSs acquire average kick velocities between 200 and 500 km s⁻¹ during their birth in supernova (SN) explosions (e.g., Helfand et al. 1977; Tutukov et al. 1984; Harrison et al. 1993; Lyne & Lorimer 1994; Kaspi et al. 1996; Hansen & Phinney 1997; Cordes & Chernoff 1998; Fryer et al. 1998; Lai et al. 2001; Arzoumanian et al. 2002; Chatterjee et al. 2005; Hobbs et al. 2005). The most likely explanation of these natal kick velocities is a recoil of the newly formed NS due to anisotropic mass ejection when the core-collapse SN explosion is initiated in a nonspherical way. The corresponding acceleration of the NS can be understood by the "gravitational tugboat mechanism," in which the asymmetrically expelled matter exerts a long-lasting momentum transfer to the NS by hydrodynamic and gravitational forces over a period of many seconds (Scheck et al. 2004, 2006; Nordhaus et al. 2010a, 2012; Wongwathanarat et al. 2010b, 2013; Bruenn et al. 2016; Janka 2017b; Müller et al. 2017).

This hydrodynamic kick scenario is supported by recent observational analyses of SN remnants with associated NSs (Bear & Soker 2017; Holland-Ashford et al. 2017; Katsuda et al. 2018), which confirm the prediction by the gravitational tugboat mechanism that the bulk of the SN ejecta should be found moving in the hemisphere opposite to the direction of the NS kick velocity (Wongwathanarat et al. 2013, 2017). This anti-alignment is imprinted in the early moments of the explosion on the chemical yields (from Si to the Fe group) that are produced in the immediate surroundings of the newborn NS, and thus are most strongly affected by anisotropic shock propagation. The asymmetry becomes more prominent with higher NS kicks. Aging SN remnants, however, will develop a global hemispheric asymmetry due to the changing relative motion of the NS and the ejecta. The SN ejecta are gradually decelerated by swept-up circumstellar material while the NS moves within the gaseous ejecta cloud with essentially maintained speed. Discriminating such an apparent kick-ejecta anti-orientation and possible effects due to environmental asymmetries from intrinsic explosion and ejecta asymmetries becomes an increasingly challenging task for older SN remnants (Katsuda et al. 2018).

Two-dimensional (axisymmetric; 2D) models (Scheck et al. 2004, 2006; Nordhaus et al. 2010a, 2012; Bruenn et al. 2016) as well as three-dimensional (3D) simulations (Wongwathanarat et al. 2010b, 2013; Müller et al. 2017) for iron-core progenitors of 15–20 M_⊙ have shown that the gravitational tugboat mechanism can explain NS kicks up to more than 1000 km s⁻¹, because considerable dipolar asymmetry in the SN ejecta can result from violent mass motions by convective overturn in the neutrino-heating layer, from large-amplitude sloshing and spiral activity of the stalled SN shock due to the standing accretion-shock instability (SASI) (e.g., Foglizzo 2002; Blondin et al. 2003; Foglizzo et al. 2006, 2007; Scheck et al. 2008) or from long-lasting, anisotropic accretion downdrafts around the nascent NS after the onset of the SN blast (Müller et al. 2017). Depending on stochastic variations of such large-scale asymmetry, the kick velocities were found to range from less than 10 km s⁻¹ to over 1000 km s⁻¹, with average values of several 100 km s⁻¹.

Apart from such stochastic case-to-case variations, hydrodynamic NS kicks associated with SN ejecta asymmetries are expected to exhibit systematic dependences on the explosion and progenitor properties. On average, larger kicks should be the result of higher explosion energies and of higher masses of the anisotropically ejected matter that accelerates the NS by momentum transfer (Janka 2017b). The kick-relevant mass involved in this process depends on the "compactness," ξ_M = (M/M_⊙)/(R(M)/1000 km), of the progenitor (defined in O'Connor & Ott (2011)), which is correlated with the steepness of the density profile around the degenerate stellar core at the onset of collapse (Sukhbold & Woosley 2014). Progenitors with shallow density profiles possess high compactness values. If they explode by the neutrino-driven mechanism, a lot of matter is available to absorb energy from neutrinos, allowing for high explosion energies and large amounts of kick-mediating ejecta. This relation explains the trend to higher NS kicks (on average) for successfully exploding high-compactness progenitors (Janka 2017b). The opposite tendency should apply for progenitors with steep density profiles and correspondingly lower values of the core compactness. Indeed, a small set of 2D and 3D explosion simulations by Suwa et al. (2015) and Müller et al. (2018) for ultrastripped SN Ic progenitors seems to support these theoretically expected relations.

In this paper, we present a large sample of 45 explosion calculations in 2D and 3D for neutrino-powered electron-capture SNe (ECSNe) of an O–Ne–Mg-core progenitor. Our results reflect the described trends much more prominently. Because of the extremely tenuous H/He envelope around the degenerate core, none of our calculations produced an NS kick of more than a few km s⁻¹. This result is in line with considerations by Podsiadlowski et al. (2004), who hypothesized that ECSNe give birth to NSs with low kick velocities. However, detailed analyses of hydrodynamic explosion models (Scheck et al. 2004, 2006; Nordhaus et al. 2010a; Wongwathanarat et al. 2010b, 2013) have revealed a different mechanism than that argued by Podsiadlowski et al. (2004). Scheck et al. (2004) found that NS kicks are caused by asymmetries in the mass ejection of SN explosions rather than anisotropies in the neutrino emission. In ECSNe, the low explosion energy and the extremely rapid outward expansion of SN shock and postshock layer (all of which are linked to the steep density decline outside of the degenerate core) are detrimental to the possibility of large NS kicks.³

Our results have important implications for the discussion of the Crab SN and its stellar origin, for which an O–Ne–Mg core progenitor has been proposed (Hillebrandt 1982; Nomoto et al. 1982). Such a connection seems to be compatible with the He-rich chemical composition of the nebula and a small amount of ejected iron-group material (Nomoto et al. 1982; MacAlpine & Satterfield 2008; Wanajo et al. 2011; Smith 2013), a small ejecta mass (Davidson & Fesen 1985; Fesen et al. 1997), low kinetic energy (Dessart et al. 2006; Kitaura et al. 2006; Fischer et al. 2010; Yang & Chevalier 2015; Radice et al. 2017), and the characteristics of the light curve (Smith 2013; Tominaga et al. 2013). However, the small NS kick velocities of only a few km s⁻¹ obtained for ECSN explosions in our work are in conflict with the ∼160 km s⁻¹ deduced from observations for the spatial velocity of the Crab pulsar (Hester 2008; Kaplan et al. 2008). This discord points to several possible consequences: If the Crab Nebula was born in an ECSN, the motion of the Crab pulsar either requires a non-hydrodynamic acceleration mechanism—for example, anisotropic neutrino emission or the post-natal electromagnetic rocket effect (Harrison & Tademaru 1975)—or the pulsar's motion is attributed to the disruption of a close binary system during the SN explosion (Blaauw effect; Blaauw 1961). In the latter case, however, the structure and internal dynamics of the Crab Nebula with its pulsar may have to be revised. Alternatively, the mismatch between our calculations for ECSNe and the observed space velocity of the Crab pulsar might imply that Crab is a relic of the explosion of a low-mass iron-core progenitor. Because such progenitors can possess considerably higher values of the core compactness (i.e., shallower density profiles around the degenerate core) than O–Ne–Mg-core progenitors, the hydrodynamical kick mechanism might account for the proper motion of the Crab pulsar in this case. We will critically assess these possibilities in a detailed discussion.

In Section 2, we summarize our numerical methods and input physics. In Section 3, we describe our set of computed models and explosion results. In Section 4 , we present our results for the corresponding NS kicks. In Section 5, we discuss the implications for Crab, and Section 6 contains a summary and our conclusions.

2.1. Numerical Code and Input Physics

2.2. Numerical Grid and Boundary Conditions

Our 2D simulations are performed with a polar coordinate grid, whereas for the 3D simulations we employ an axis-free Yin-Yang grid (Kageyama & Sato 2004), implemented into Prometheus-HotB by Wongwathanarat et al. (2010a) as an overlay of equatorial patches of two spherical polar grids, tilted by 90 degrees relative to each other. The Yin-Yang grid avoids the coordinate singularity along the polar axis of the spherical grid, and thus permits computations with considerably higher efficiency by relaxing the Courant–Friedrichs–Lewy (CFL) condition on the time step.

We use 1310 zones in the radial direction for both 2D and 3D models. For all simulations, we employ the same non-equidistant spacing of the radial grid points, similar to the grid described by Janka et al. (2008), in order to resolve the steep density gradient at the surface of the O–Ne–Mg core. In the angular directions, the grid is chosen to be equidistant in both 2D and 3D. We apply an angular resolution of 1° for the 2D simulations and of 3° for our 3D simulations. This corresponds to 180 azimuthal zones in 2D and 30(θ) × 90(ϕ) angular zones on each of the two Yin-Yang patches of the 3D grid.

In both the 2D and 3D calculations, we exclude the high-density, neutrino-opaque core of the newly formed NS and replace it by a gravitating point mass at the grid center, as well as an inner, moving grid boundary at a prescribed, time-dependent radius, following the numerical treatment used already in previous publications (Janka & Müller 1996; Kifonidis et al. 2003; Scheck et al. 2006; Arcones et al. 2007). Thus excising the central region of the proto-NS (PNS) from the hydrodynamical domain relaxes the time-step constraint. Therefore, it permits us to compute a large set of 2D and 3D models for one second of physical evolution time with reasonable computational resources. A full treatment from first principles would have been computationally prohibitive, particularly in 3D, for gathering model statistics over a larger sample of simulation runs. Moreover, the grid boundary provides us with a simple handle to regulate the SN explosion energy to preferred numbers by choosing suitable values of the free parameters involved in our description of the boundary condition.

Such a (relatively) free choice of the explosion energy⁴ in our study is desirable not only because it allows us to account for uncertainties of the ECSN energy on the observational side (e.g., Tominaga et al. 2013; Yang & Chevalier 2015), but also for remaining uncertainties in the first-principle modeling of ECSN explosions. Theoretical uncertainties are, for example, connected to the still incompletely determined properties of the high-density EoS in hot NSs and the corresponding neutrino opacities. Both have an influence on the exact energetics of the neutrino-driven wind of the PNS, and hence also on the energy of the neutrino-powered SN blast (von Groote 2014). For our purpose, the use of a parameterized model to avoid a full numerical treatment of the high-density NS core is therefore preferable because we are interested in a large, representative study of the multidimensional dynamics of the SN ejecta for defined explosion energies. Quantities that are directly influenced by the physics within the nascent NS, such as the exact properties of the emitted neutrino signal or the internal evolution of the NS itself, are not of interest for our study.

Following Scheck et al. (2006) and Arcones et al. (2007), and motivated by the shrinking of the NS core as it cools and deleptonizes by neutrino emission, we apply a contraction of the inner-boundary radius:

$$\begin{eqnarray}&&{R}_{\mathrm{ib}}(t)={R}_{\mathrm{ib}}^{{\rm{f}}}+({R}_{\mathrm{ib}}^{{\rm{i}}}-{R}_{\mathrm{ib}}^{{\rm{f}}})\exp (-t/{t}_{\mathrm{ib}}).\end{eqnarray} \tag{ 1 }$$

This prescription causes the inner grid boundary to retreat from an initial radius ${R}_{\mathrm{ib}}^{{\rm{i}}}$ to a final radius ${R}_{\mathrm{ib}}^{{\rm{f}}}$ with a characteristic exponential timescale t_ib. For the progenitor studied in this work, we choose ${R}_{\mathrm{ib}}^{{\rm{i}}}=43.5$ km. This corresponds to an enclosed baryonic mass of ${M}_{\mathrm{encl}}^{{\rm{b}}}=1.02$ M_⊙ in the post-bounce initial state from which we start our multidimensional simulations. This enclosed baryonic mass is kept constant with time, for which reason R_ib(t) tracks a fixed mass coordinate. The hydrodynamic quantities that must be defined in the ghost cells at the inner boundary are chosen such that the condition of hydrostatic equilibrium is fulfilled (Janka & Müller 1996).

At the outer boundary of the computational mesh, at a radius of R_ob = 4.6 × 10⁵ km, a free in/outflow boundary condition is applied. This radius is picked such that the SN shock does not leave the computational domain within the simulated evolution.

As in Scheck et al. (2006) and Arcones et al. (2007), we impose spherically symmetric luminosities for all neutrino species at the inner grid boundary.⁵ These core luminosities are regulated by a parameter ${\rm{\Delta }}{E}_{\nu ,\mathrm{core}}^{\mathrm{tot}}$, which determines the total gravitational binding energy extracted by neutrino radiation from the NS core. In the present study, we use three parameters of the NS-core model, namely ${\rm{\Delta }}{E}_{\nu ,\mathrm{core}}^{\mathrm{tot}}$, ${R}_{\mathrm{ib}}^{{\rm{f}}}$, and t_ib (Table 1), to regulate the explosion energies of our computed ECSN models (see Section 3). Although these parameters have the physical meaning mentioned above and discussed in detail by Scheck et al. (2006), their values are varied here beyond the ranges suggested by self-consistent simulations of the NS-core evolution. We therefore discourage a deeper interpretation of the numbers given in Table 1; they are picked merely with the goal of producing desired values of the SN explosion energy by neutrino heating. The corresponding combinations of parameter values are not unambiguous, i.e., other choices can be equally suitable to obtain the same explosion energies.

Table 1.  Parameters Values for the Five Model Sets with Chosen SN Energies

Model Set	Eexp	${\rm{\Delta }}{E}_{\nu ,\mathrm{core}}^{\mathrm{tot}}$	${R}_{\mathrm{ib}}^{{\rm{f}}}$	tib
	(1049 erg)	(1051 erg)	(km)	(s)
O3	∼3	0.2	20	2.0
O5	∼5	3.6	15	1.0
O8	∼8	17.9	15	0.7
O12	∼12	35.8	15	0.5
O16	∼16	53.6	17	0.3

Download table as:  ASCIITypeset image

2.3. Seed Perturbations

3.1. Model Sets

We investigate the explosion dynamics of an 8.8 ${M}_{\odot }$ star with an O–Ne–Mg core prior to collapse (Nomoto 1984, 1987); plots of the progenitor structure, also in comparison to iron-core models and more modern O–Ne–Mg core progenitors, can be found in Janka et al. (2008, 2012), Jones et al. (2013), Müller (2016), Janka (2017a), and Cerdá-Durán & Elias-Rosa 2018).⁶ Rotation is not accounted for in our study.⁷

The collapse and early post-bounce evolution of the considered ECSN progenitor were simulated in spherical symmetry by Hüdepohl et al. (2010) using the Prometheus-Vertex code (Rampp & Janka 2002), which, unlike Prometheus-HotB, includes the whole core without parameterized inner boundary condition. It also contains all relevant microphysics (in particular, the electron-capture rates and a special treatment of highly degenerate conditions) that are needed to follow the core-collapse and shock-formation phases accurately. At 18 ms post-bounce, we adopt the hydrodynamic model state from the Prometheus-Vertex calculation and track the subsequent evolution via multidimensional simulations with the Prometheus-HotB code. The mapping from one simulation to the other, including a change of the numerical grid (which has essentially the same radial resolution in both cases) and of the neutrino treatment does not produce any noticeable numerical artifacts. This allows us to perform 40 explosion runs in 2D, and five more models in 3D, until 1 s physical time. Note that, in the whole of our paper, unless stated otherwise, time is always defined by the simulated time of evolution after taking over the computation with Prometheus-HotB. For the true post-bounce time, 18 ms have to be added to this clock. The large number of simulations enables us to explore the statistical fluctuations of the NS kicks.

The free parameters that describe the core-neutrino source in our setup, ${R}_{\mathrm{ib}}^{{\rm{f}}}$, t_ib, and ${\rm{\Delta }}{E}_{\nu ,\mathrm{core}}^{\mathrm{tot}}$ (see Section 2.2), are calibrated such that our models reproduce basic explosion features of ECSNe found in fully self-consistent simulations (Kitaura et al. 2006; Janka et al. 2008; Fischer et al. 2010; Hüdepohl et al. 2010; von Groote 2014; Radice et al. 2017), particularly the explosion energy and the time of the onset of the SN blast. We choose five parameter triples for covering the possible range of explosion energies suggested by the mentioned full-scale ECSN models and compatible with SN 1054 according to a detailed analysis of the Crab remnant (Yang & Chevalier 2015). In Table 1, the corresponding values of the three relevant calibration parameters are listed for our five model sets.

We use the following naming convention for the computed models in these five sets (see also Table 2): The initial letter "O" stands for "O–Ne–Mg core" and emphasizes that we consider the same progenitor star for all calculations. The associated number from the set {3, 5, 8, 12, 16} represents the approximate value of the explosion energy (in units of 10⁴⁹ erg). For the 2D models, the next two letters specify how the initial model is randomly perturbed (cell-by-cell) at the beginning of our simulation in order to break spherical symmetry and thus to initiate the growth of nonradial hydrodynamic instabilities. The first letter represents the perturbed quantity, i.e., either the radial velocity v_r ("v") or the density ρ ("d"). In the case of radial velocity perturbations, conservation of the total (i.e., internal plus kinetic) energy is enforced. In the case of density perturbations, conservation of the total mass is ensured. The amplitude applied for the initial perturbation is indicated by the second letter: "h" (high), "i" (intermediate), and "l" (low) for an amplitude of 1%, 0.5%, or 0.1%, respectively. Finally, models computed for identical parameter settings are numbered sequentially. To facilitate orientation, we summarize our naming convention in Table 2. For all of the models computed in 3D, we apply an initial perturbation of amplitude 0.1% in the radial velocity. These models can be identified by the suffix "3D" in the model name. In our 2D simulations, both velocity and density profiles are considered for perturbations with various amplitudes, in order to verify that the evolving stochasticity of the hydrodynamic instabilities is basically robust to different types of perturbations. Indeed, we do not observe any qualitative differences in the outcomes caused by different quantities and patterns of perturbation. For our 3D models, we cannot provide the same variety of cases as for the 2D models, and employ merely one perturbation setting for all 3D simulations. A "set of models" is defined by all explosion simulations that share the same values of the calibration parameters, and thus develop very similar explosion energies, independent of their specific perturbation pattern.

Each of our five model sets consists of five simulation runs. In other words, for each combination of values for our triple of parameters, we performed five simulations using different seed perturbations. One corresponding case per set was computed in 3D. Model set O5 is treated preferentially by repeating the 2D calculations for these parameter values 20 times. The reason for this exception is that a recent detailed analysis of the Crab Nebula tends to favor a low blast-wave energy (around 5 × 10⁴⁹ erg) for SN 1054 (Yang & Chevalier 2015), in line with fully relativistic 1D and 2D calculations of O–Ne–Mg core explosions (von Groote 2014).

All simulations were carried on for 1 s of physical evolution, except for O16-3D, which could not be advanced beyond 750 ms. At that time, this run terminated because of insufficient resolution at the shock at very large distances from the center. Because this problem does not affect the NS-kick relevant inner ejecta and all characteristic parameters of the explosion (in particular, blast-wave energy, NS velocity, and acceleration) are basically converged, we accept the output at 750 ms as the final result for model O16-3D.

Figure 1 displays the time evolution of the diagnostic explosion energies for all simulated 2D and 3D models. The variability for a fixed parameter set is caused by stochastic effects in the hydrodynamic evolution as a result of the nonlinear growth of hydrodynamic instabilities triggered by the initial random perturbations. Moreover, Figure 1 introduces a color coding for the model of our five sets that will also be applied in Figures 8 and 9. The diagnostic energy is defined as the volume integral over the sum of internal, kinetic, and gravitational energy densities for all post-shock grid cells where this sum is positive.⁸ Shortly after the SN shock has begun to accelerate outward, the latter condition is fulfilled for all expanding, high-entropy material in the layer between shock front and gain radius, i.e., in the region where neutrino heating is responsible for net energy deposition. Once the shock has reached a radius of several 1000 km, the diagnostic energy is essentially equal to the (instantaneous) explosion energy, because the gravitational binding energy of the overlying stellar layers ahead of the shock becomes extremely low in the considered progenitor (see Figure 1 in Radice et al. 2017). Because the final explosion energy is basically reached after ∼300 ms at the latest, we will not distinguish between diagnostic energy and explosion energy in the remainder of this paper.

Zoom In Zoom Out Reset image size

3.2. Explosion Properties of ECSNe

Our simulations exhibit an explosion behavior very similar to the published ECSN models discussed in the literature: see Kitaura et al. (2006), Fischer et al. (2010), and Hüdepohl et al. (2010) for 1D simulations, and see Janka et al. (2008), Wanajo et al. (2011), von Groote (2014), and Radice et al. (2017) for 2D results. We therefore summarize here only the basic features and focus specifically on those aspects that are of interest in the context of the NS kicks that will be discussed in the subsequent sections.

The extremely steep density gradient at the edge of the O–Ne–Mg core crucially determines the shock dynamics in ECNSe. In response to the rapidly decreasing ram pressure of the infalling matter when the shock reaches the core/envelope boundary, fast outward shock acceleration sets in after less than ∼100 ms of slower—but continuous—shock expansion (very similar to the evolution shown in Figure 3 of Janka et al. 2008). Within this time, nonradial asymmetries begin to develop only in the close vicinity of the NS, because neutrino-heated high-entropy matter outside of the gain radius becomes buoyant and starts rising in Rayleigh–Taylor plumes. This phenomenon takes place basically independent of the explosion energy and of the dimension of the models (see the upper row of Figure 2, and the upper left panel of Figure 3). When the most extended of these convective bubbles reach radii of a few hundred kilometers at ∼100 ms, the shock has already propagated to more than 10⁴ km (see Figure 3 in Janka et al. 2008) and is not affected by the hydrodynamic instabilities growing far behind in its wake. An almost spherical expansion of the shock is therefore a generic feature of ECSNe.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 2 shows snapshots of the 2D entropy distribution at different times between 100 ms and 1 s for the axisymmetric models O5-vl (left column) and O16-vl (right column), which represent a low-energy and a high-energy explosion, respectively. Both models exhibit an almost self-similar expansion of the Rayleigh–Taylor plumes that have developed already until 100 ms. The expansion is slightly faster for the explosion of model O16-vl, which is over three times more energetic.

The lower-energy model O5-vl and the higher-energy model O16-vl are representatives of two antipodes with respect to their behaviors. In O16-vl, the neutrino luminosities from the PNS are high enough to drive a powerful baryonic outflow ("neutrino wind") via neutrino-energy deposition near the PNS surface. This essentially spherically symmetric mass outflow leads to the continuous, slow growth of the explosion energy visible in Figure 1 from ∼400 ms until the end of the simulation. During this phase, the high-entropy wind region can be seen as a red sphere around the center in the right panels of Figure 2. Because the wind fills the central volume, the Rayleigh–Taylor plumes of the convective layer are pushed outward and exhibit signs of a growing compression in the radial direction. The two distinct downflows at the poles (with the stronger one at the north pole along the +z-direction) are axis artifacts fostered by the flow constraints connected to the assumption of axisymmetry.

In contrast, the neutrino heating by the lower luminosities in model O5-vl is too weak to shed matter from the PNS surface in a strong neutrino-driven wind. Therefore, the explosion energy levels off to its final value after ∼300 ms (Figure 1), and convective filaments keep billowing around the NS until the end of the simulation. The initial Rayleigh–Taylor plumes continue to expand in a self-similar manner.

Figures 3 and 4 depict the entropy per nucleon, color coded, at different times in a cross-sectional plane of the 3D data from the evolution of model O8-3D. This model develops the largest ejecta asymmetry and the highest NS kick of all our 3D runs. The angular scale of the Rayleigh–Taylor plumes is corroboratively similar to the 2D cases displayed in Figure 2 and clearly dominated by higher-order spherical harmonics and nearly self-similar expansion already after 150 ms (see Figure 5). The sequence of plots in Figure 4 visualizes the formation of a prominent, long-lasting downflow, which grows after ∼400 ms from a massive pocket of cool, low-entropy material enclosed by the outward-rising bubbles. While the development of buoyant Rayleigh–Taylor plumes of neutrino-heated matter with similar geometrical patterns is a generic feature of all of our 3D models, the long-lasting downflow is specific to model O8-3D. This downflow determines the final direction and magnitude of the NS kick, which is generated by the gravitational tugboat mechanism and is largest in model O8-3D, compared to the other 3D cases, to be discussed in the next section.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.1. NS Kicks—Physics

The gravitational tugboat mechanism attributes natal kicks of NSs to anisotropic ejection of matter that exerts gravitational and hydrodynamic forces on the NS (Scheck et al. 2006; Wongwathanarat et al. 2010b, 2013; Janka 2013, 2017b). In our simulations, the NS is tied to the center of the numerical grid because of the excluded core volume, which is replaced by a gravitating point mass and an inner grid boundary. However, due to the asymmetries growing via nonradial hydrodynamic instabilities at the beginning of the SN explosion, the ejecta can acquire a net linear momentum on the grid. The corresponding recoil velocity of the NS, ${{\boldsymbol{v}}}_{\mathrm{ns}}$, can be evaluated by applying linear momentum conservation of the ejected gas plus the compact remnant, i.e., the momentum carried by the NS is given by the negative of the ejecta momentum ${{\boldsymbol{P}}}_{\mathrm{gas}}$. This yields

$$\begin{eqnarray}&&{{\boldsymbol{v}}}_{\mathrm{ns}}(t)=-\displaystyle \frac{{{\boldsymbol{P}}}_{\mathrm{gas}}(t)}{{M}_{\mathrm{ns}}(t)},\end{eqnarray} \tag{ 2 }$$

where the mass of the NS is given by M_ns, which we take to be the baryonic mass enclosed by the NS radius R_ns. The latter is defined as the (angle-dependent) radius where the density drops below 10¹¹ g cm⁻³. The gas momentum ${{\boldsymbol{P}}}_{\mathrm{gas}}={\int }_{{R}_{\mathrm{ns}}}^{{R}_{\mathrm{ob}}}\rho {\boldsymbol{v}}\,{dV}$ is the net linear momentum on the computational grid exterior to the NS. It should be kept in mind that this momentum is exclusively associated with the ejecta mass, ${M}_{\mathrm{ej}}={\int }_{{R}_{\mathrm{ns}}}^{{R}_{{\rm{s}}}}\rho \,{dV}$ (with R_s being the angle-dependent shock radius), because the linear momentum of the gas ahead of the shock is negligible (given that the shock is nearly spherical around the grid center and the star is spherical and basically at rest).

In order to better understand the exact acceleration mechanism acting on the NS and to support the direct numerical evaluation of the simulation results with a more detailed physical consideration (following previous works regarding NS kicks in SNe from Fe-core progenitors, e.g., Scheck et al. (2006), Nordhaus et al. (2010a), Wongwathanarat et al. (2010b, 2013), Müller et al. (2017)), we separate the individual forces that act on the NS and cause its kick (Scheck et al. 2006):

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d}{{dt}}{{\boldsymbol{P}}}_{\mathrm{ns}} & = & \displaystyle \frac{d}{{dt}}({M}_{\mathrm{ns}}\,{{\boldsymbol{v}}}_{\mathrm{ns}})\\ & = & -{\oint }_{r={r}_{0}}p\,d{\boldsymbol{S}}-{\oint }_{r={r}_{0}}\rho {\boldsymbol{v}}{v}_{r}\,{dS}+{\displaystyle \int }_{r\gt {r}_{0}}{{GM}}_{\mathrm{ns}}\displaystyle \frac{{\boldsymbol{r}}}{{r}^{3}}\,{dm}\end{array}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\approx {M}_{\mathrm{ns}}\,\displaystyle \frac{d{{\boldsymbol{v}}}_{\mathrm{ns}}}{{dt}},\end{eqnarray} \tag{ 4 }$$

where ${dS}={r}^{2}\sin \theta d\theta d\phi$ is the surface element on a sphere of radius r. Using the approximative relation in the third line of this equation (and thus adopting the analysis of NS kicks in Fe-core SNe), dividing by M_ns, and integrating over time yields the NS recoil velocity, provided M_ns does not change with time. For the calculation of NS kicks in SNe of Fe-core progenitors, this turned out to work very well, because in those cases the crucial phase of NS acceleration occurs after the onset of the SN explosion and continues for several seconds; see Scheck et al. (2006), Wongwathanarat et al. (2010b, 2013), Nordhaus et al. (2010a), and Müller et al. (2017). During this phase of evolution, the assumption of a constant PNS mass is good roughly on the percent level. In the case of the O–Ne–Mg core explosions considered here, however, the NS kicks happen so fast and early that variations of the NS mass might have a more relevant impact. We will discuss this in Section 4.2.

The first term on the r.h.s. of Equation (3) accounts for pressure variations over a spherical surface chosen at r₀ = 1.3R_ns (Wongwathanarat et al. 2013). The momentum flux across this surface (second term) can be decomposed into down- and outflows, depending on the sign of the radial velocity v_r. Both terms only contribute significantly if anisotropically distributed matter is in sonic contact with the sphere of evaluation. Otherwise, if asymmetries are present only in matter that is hydrodynamically detached from the NS, the only remaining possibility to accelerate the NS is via the long-range forces of gravity, which are accounted for by the third term. In contrast to the hydrodynamic forces, the NS interacts gravitationally with matter in the entire volume under consideration. Therefore, the third term typically constitutes the longest-lasting contribution to the total acceleration. As the gravitational force drops with the square of the distance, material that stays longer in the vicinity of the NS has a larger accelerating influence on it. This so-called gravitational tugboat mechanism (Wongwathanarat et al. 2013) pulls the NS into the direction of the most massive, slowest-moving ejecta. These ejecta preferentially lie in the direction where the explosion is weaker. This means that the NS is kicked opposite to the hemisphere of the stronger explosion, as expected for a hydrodynamical acceleration associated with blast-wave asymmetries.

4.2. Results

4.2.1. Exemplary Cases

With insight gained about the forces that contribute to the acceleration of the NS, we return to the exemplary 2D models O5-vl and O16-vl shown in Figure 2. In such calculations, the assumed axisymmetry constrains the NS kick direction to coincide with the symmetry axis. For both models, we obtain a positive kick velocity, i.e., a kick vector pointing to the (+z)-direction (see Table 3 for results of all models). This complies with the observation that both models develop a stronger downflow at the north pole, indicating the hemisphere where the SN blast is weaker. While O16-vl is on the extreme side, with its NS kick velocity of more than 6 km s⁻¹, model O5-vl is a fairly average 2D case, with ∼3 km s⁻¹ at the end of our simulation.

Table 3.  Final Values of Explosion and NS Properties for All 2D and 3D Models^a

Model	Eexp	texp	Mns	$\langle {v}_{\mathrm{ns}}\rangle$	$\langle {a}_{\mathrm{ns}}\rangle$	Mej	Pej	$\langle {\alpha }_{\mathrm{ej}}\rangle$
	(1049 erg)	(ms)	(${M}_{\odot }$)	($\mathrm{km}\,{{\rm{s}}}^{-1}$)	(km s−2)	(10−2 ${M}_{\odot }$)	(1040 cm g s−1)	(%)
O3-dl	3.60	106	1.364	−1.4	1.6	1.24	3.26	1.2
O3-vh	3.78	102	1.364	−0.1	−0.4	1.23	3.31	0.1
O3-vi	3.68	102	1.364	−2.2	−0.4	1.25	3.26	1.8
O3-vl	2.69	132	1.365	−0.0	−0.3	1.17	2.68	0.0
O3-vl2	2.73	130	1.365	0.8	−1.2	1.15	2.68	0.8
O5-dl	4.86	110	1.363	−0.1	0.2	1.33	4.35	0.1
O5-vh	5.70	96	1.363	2.1	0.7	1.32	4.69	1.2
O5-vi	5.90	98	1.363	5.2	1.2	1.39	4.91	2.9
O5-vl	4.84	110	1.364	2.9	1.5	1.26	4.18	1.9
O5-vl2	5.12	106	1.364	2.7	−0.5	1.25	4.31	1.7
O5-vl3	5.13	116	1.363	−0.1	−2.0	1.36	4.52	0.1
O5-vl4	5.23	108	1.364	5.8	0.5	1.27	4.38	3.6
O5-vl5	4.85	110	1.364	3.2	−1.5	1.27	4.19	2.0
O5-vl6	4.96	108	1.364	−0.7	0.4	1.23	4.19	0.5
O5-vl7	5.18	106	1.364	−0.7	−1.4	1.31	4.45	0.4
O5-vi2	5.67	98	1.363	1.3	−1.1	1.39	4.78	0.7
O5-vi3	5.90	98	1.362	1.0	0.2	1.44	4.99	0.5
O5-vi4	5.79	98	1.363	0.8	0.8	1.36	4.80	0.5
O5-vi5	5.86	100	1.363	0.4	1.7	1.41	4.89	0.2
O5-vi6	6.03	98	1.362	−1.9	−0.1	1.46	5.07	1.0
O5-vh2	5.71	98	1.363	0.4	0.8	1.33	4.74	0.2
O5-vh3	6.53	98	1.361	0.3	0.2	1.54	5.49	0.1
O5-vh4	5.49	98	1.363	0.9	0.9	1.33	4.61	0.5
O5-vh5	6.17	96	1.362	1.2	0.1	1.44	5.09	0.7
O5-vh6	6.07	98	1.363	1.3	0.3	1.40	5.00	0.7
O8-dl	9.01	98	1.359	−4.9	−1.0	1.75	7.09	1.9
O8-vh	8.96	94	1.359	−2.0	−0.8	1.77	7.07	0.8
O8-vi	8.89	96	1.359	4.0	1.6	1.74	6.97	1.5
O8-vl	8.36	98	1.360	−0.6	0.2	1.72	6.69	0.3
O8-vl2	8.52	100	1.360	2.4	1.1	1.70	6.77	1.0
O12-dl	12.16	96	1.356	−2.1	−0.2	2.05	9.05	0.6
O12-vh	12.64	92	1.356	−3.7	−0.3	2.07	9.25	1.1
O12-vi	12.38	94	1.356	−1.5	−0.8	2.02	9.10	0.4
O12-vl	12.25	96	1.356	−0.4	−0.1	2.05	9.09	0.1
O12-vl2	12.23	96	1.356	−2.5	−0.2	2.04	9.04	0.7
O16-dl	16.55	92	1.353	2.7	0.1	2.41	11.53	0.6
O16-vh	16.85	90	1.352	0.6	0.0	2.41	11.59	0.1
O16-vi	16.92	88	1.353	−3.0	−0.2	2.37	11.51	0.7
O16-vl	16.42	92	1.353	6.3	0.2	2.40	11.45	1.5
O16-vl2	16.59	92	1.353	−3.4	−0.1	2.39	11.48	0.8
O3-3D	3.41	111	1.364	0.5	0.0	1.24	3.19	0.4
O5-3D	5.73	105	1.362	0.9	1.4	1.48	4.91	0.5
O8-3D	8.35	108	1.359	2.5	1.2	1.75	6.74	1.0
O12-3D	12.00	107	1.356	0.4	−0.4	2.08	9.04	0.1
O16-3Db	15.92	101	1.353	0.3	−0.3	2.38	11.04	0.1

Notes.

^aSee Section 4.2.3 for definitions. ^bSimulation terminated after 0.75 s physical time, all other models were evolved for 1 s.

Download table as:  ASCIITypeset image

Because of the development of one distinct downflow during the later stages (t ≳ 400 ms) of the computed evolution, our example 3D model O8-3D (Figures 3 and 4) exhibits a kick geometry very similar to the discussed 2D cases. The kick direction almost coincides with the downflow, indicating that this downflow is the main feature that determines the NS kick direction and magnitude. The deviation from a perfect alignment is caused by the massive ejecta clumps that expand with relatively lower velocities in the right hemisphere (yellow regions with intermediate entropy values and relatively high densities at the one o'clock and five o'clock positions in the bottom right panel of Figure 4). Model O8-3D, with its prominent single downflow, however, is not representative of our sample of computed 3D cases. Its NS kick velocity of 2.5 km s⁻¹, which is still growing at a comparatively high rate of 1.2 km s⁻² at one second, is the largest of all our 3D models. This suggests that the ejecta geometry that has formed in model O8-3D is not a generic feature in 3D, as mentioned before (Section 3.2). The relation of explosion asymmetries and NS kicks in 2D and 3D models with be further discussed in Sections 4.2.2 and 4.2.3.

Figures 6 and 7 display the time evolution of the NS velocities and the different contributions to the NS acceleration computed with Equations (2)–(4) for the models under consideration. In Figure 6, positive values correspond to motion/acceleration in the (+z)-direction of the 2D polar grid. The NS velocities computed directly from the hydrodynamical results with Equation (2), under the assumption of momentum conservation (upper panels; black lines), agree very well with the summed effects of all accelerating forces as obtained by an analysis using Equation (3) (upper panels; red lines). The less energetic explosion of model O5-vl features stronger accelerations, and the NS velocity exhibits more dramatic changes between positive and negative values during the initial ∼200 ms than in the higher-energy explosion of model O16-vl, where the NS velocity is fairly monotonically boosted to its final value. In both cases, the terminal speed of the NS is reached within roughly 300 ms. However, in model O5-vl, the NS motion continues to exhibit low-amplitude fluctuations caused by the persisting nonradial flows in the low-density medium around the NS (see Figure 2, left column). These lead to variations of the acceleration, the amplitudes of which decrease only slowly with time.

Zoom In Zoom Out Reset image size

In model O8-3D (Figure 7), the magnitude and direction of ${{\boldsymbol{v}}}_{\mathrm{ns}}$ have not yet asymptoted by the end of the simulation. However, the NS acceleration decreases continuously, so a tendency toward saturation is clearly visible and a further growth of the NS kick by at most ∼0.5–1 km s⁻¹ may be expected. The velocity components in the (−y) and (−z)-directions dominate, as may be concluded from an inspection of the orientation of the kick vector relative to the tripod in Figure 4.

Zoom In Zoom Out Reset image size

The dominant contributions to the NS acceleration result from gravitational and pressure forces (orange and blue lines, respectively, in the acceleration panels of Figures 6 and 7), which essentially always counteract and nearly balance. Because pressure forces require matter to be in sonic contact with the NS surface, the corresponding acceleration is exerted by mass in the very close vicinity of the NS. Gravity causes almost the same effect as pressure forces, but with opposite sign, which suggests that the acceleration by gravity forces is also mostly associated with matter in the near surroundings of the NS in our models. That is, the same clumps of matter that pull the NS gravitationally also push it away by means of pressure. Because of the extremely rapid outward expansion of the ejecta in ECSNe and the lack of a large-scale dipolar asymmetry, the gravitational influence of outflowing gas quickly abates when this material rushes to larger radial distances. This explains why also the gravitational influence on the NS is dominantly linked to nearby matter.

Hydrodynamic mass flows toward as well as away from the NS always play only an ancillary role for the acceleration of the NS. The medium surrounding the NS has a fairly low density, and powerful accretion flows with supersonic velocities (like in explosions of more massive progenitors) are absent. For this reason, the momentum fluxes associated with anisotropic flows around the NS are very small. The very fast expansion of the asymmetric ejecta and the rapid decline of the densities around the PNS in ECSNe also explain why the NS acceleration is finished within only ∼300 ms, i.e., a much shorter timescale than in massive Fe-core progenitors, where the NS experiences an acceleration over several seconds (Scheck et al. 2006; Wongwathanarat et al. 2013; Müller et al. 2017).

Note that minor differences between the NS velocities determined with Equations (2) and (3) in the cases of models O5-vl and O8-3D (black and and red lines in the panels of Figures 6 and 7, showing v_ns versus t) develop at an early stage of the evolution (100–300 ms p.b.) and stay basically constant later on. These differences can be understood through two facts. First, during the early post-bounce phase, mass inflows and outflows as well as gravitational settling can still lead to small changes of the PNS mass.⁹ However, when computing ${{\boldsymbol{v}}}_{\mathrm{ns}}(t)$ by integrating Equation (3), we take into account only the first term in the expression $d{{\boldsymbol{P}}}_{\mathrm{ns}}/{dt}={M}_{\mathrm{ns}}d{{\boldsymbol{v}}}_{\mathrm{ns}}/{dt}\,+{{\boldsymbol{v}}}_{\mathrm{ns}}{{dM}}_{\mathrm{ns}}/{dt}$ but ignore the second term (which leads to the approximation of Equation (4)). This implies an error of order ΔM_ns/M_ns in the corresponding estimate of the NS velocity, when ΔM_ns is the change of the NS mass during the phase of NS acceleration. In our O–Ne–Mg core explosions, M_ns varies by less than ∼5% at t ≳ 100 ms post bounce and by less than ∼1% at t ≳ 200 ms post bounce. For this reason, the omission of the second term of $d{{\boldsymbol{P}}}_{\mathrm{ns}}/{dt}$ in Equation (4) in combination with Equation (3) can cause errors only of order ∼1%–5% in $| {{\boldsymbol{v}}}_{\mathrm{ns}}|$, and even this is only at t ≲ 200 ms after bounce. The main reason for discrepancies larger than that between results from Equations (2) and (3), in particular at t ≳ 200 ms after bounce, are numerical integration errors associated with the discrete time sampling of the simulation outputs that are subsequently post-processed for evaluating the different terms in Equation (3). When the NS kick develops in a phase of very rapid variations, the finite time sampling, in steps of typically about 0.5–1 ms, can be too coarse to capture all spikes in the acceleration terms. Usually, the associated errors are of minor relevance when the net NS kicks are large or temporal variations of the acceleration terms happen on relatively long timescales. In the present models, however, the different contributions to the NS acceleration fluctuate rapidly and mostly cancel each other, particularly when the net kicks become very small. This is the case, e.g., in the left panels of Figure 6 and in Figure 7, which display cases with relatively large deviations between results from momentum conservation in the hydrodynamic simulation (Equation (2); black lines) and results from post-processing of the accelerating forces (Equation (3); red lines). Despite asymptotic quantitative differences on the level of 10%–20% in these two cases, however, both the black and red lines still follow each other closely and even match quantitatively in the right panels of Figure 6. We therefore deem the sum of the force terms considered in Equation (3) to be a satisfactory confirmation for the NS kicks deduced from our hydrodynamic simulations.

A second note at this place concerns a possible neutrino-induced kick. Evaluating our models for anisotropic neutrino emission, we obtain only tiny effects, because the lack of high-density layers around the degenerate progenitor core prevents the occurrence of massive accretion flows to the NS, at least according to our explosion models (and understanding of the explosion mechanism) and for current ECSN progenitors. Even in massive iron-core progenitors, neutrino-induced kicks due to anisotropic neutrino emission from asymmetric accretion flows are usually very small compared to the hydrodynamical kicks mediated by the gravitational tugboat mechanism. Corresponding results and the reasons for this finding are presented and explained in Scheck et al. (2006) and Wongwathanarat et al. (2013). However, a possible NS kick by anisotropic transport of neutrinos out of the NS interior cannot be evaluated in our simulations, because the NS core is excised from the computational domain and the NS evolution is not treated in detail. Instead, spherically symmetric neutrino luminosities are imposed at the contracting inner grid boundary, as described in Section 2.2. This prescription torpedoes any finally conclusive analysis of neutrino-induced NS kicks on the basis of the present models. We will return to this question again in the discussion of Section 5.

4.2.2. Compilation of Simulation Results

We evaluate the NS kick for all models according to Equation (2). The corresponding results are listed in Table 3, and Figure 8 depicts a histogram, normalized with regard to the total number of simulation runs, for the absolute values of the NS kick velocities of all computed models, binned in intervals of 0.5 km s⁻¹. None of the computed cases yields a kick larger than 7 km s⁻¹ after one second of simulated evolution, and the NS kick is less than 1.5 km s⁻¹ in more than half of the model runs. In particular, all of the 3D models tend to lie on the low-velocity side of the distribution, and the highest kick in 3D (model O8-3D) is less than half the maximum value obtained for the 2D models.

Zoom In Zoom Out Reset image size

It is clear from these results that ECSNe are unlikely to produce hydrodynamical NS kicks of more than a few kilometers per second. This is in stark contrast to explosions of massive Fe-core progenitors, where explosion models yield average kick velocities of several hundreds of km s⁻¹ by the same acceleration mechanism (Scheck et al. 2004, 2006; Nordhaus et al. 2010a, 2012; Wongwathanarat et al. 2010b, 2013; Bruenn et al. 2016; Müller et al. 2017), in agreement with observations of many radio pulsars as well as NSs in SN remnants. The NS kicks we obtain for ECSNe are even an order of magnitude lower than the kick velocities (up to some 10 km s⁻¹) found for explosions of ultrastripped SN progenitors by Suwa et al. (2015) and Müller et al. (2018).

From our discussion of the explosion dynamics of ECSNe (Section 3.2) and of the basic physics of the gravitational tugboat mechanism providing the NS acceleration (Section 4.1), we can understand the reasons for the inefficiency of ECSNe in kicking the NS. Because of the extremely steep density decline at the edge of the degenerate core, the SN shock runs outward with very high speed, also allowing the rapid expansion of neutrino-heated matter in its wake. Although convective overturn develops above the gain radius, the high-entropy plasma expands so quickly that the buoyant Rayleigh–Taylor plumes freeze out in a short-wavelength pattern. This implies that a sizable dipole mode is absent and the momentum asymmetry therefore remains small. Moreover, efficient momentum transfer from the anisotropic ejecta to the NS is prevented by the fast outward acceleration and by extremely low densities in the surroundings of the PNS, entailing that very little mass carries the asymmetries at the beginning of the explosion.

In SN explosions of massive iron-core progenitors, significant dipolar asymmetry of the neutrino-heated ejecta is a generic feature (e.g., Wongwathanarat et al. 2013; Lentz et al. 2015; Melson et al. 2015a; Müller et al. 2017; Ott et al. 2018). It originates from the development of low-mode mass motions associated with convective or SASI activity in the volume between PNS and stalled SN shock. These instability modes take about 100 ms to reach the nonlinear stage, and SASI growth is particularly fast during phases of shock stagnation or even retraction. In ECSNe, however, the extremely steep density drop outside of the O–Ne–Mg core (or equivalently, the extremely small compactness value of the core¹⁰ ) permits continuous shock expansion (even during the first 100 ms after bounce) without any phase of shock stagnation (see Janka et al. 2008). This disfavors the growth of dipolar asymmetry modes in the neutrino-heated layer as visible in Figure 5, where the Rayleigh–Taylor plumes reflect a high-order spherical harmonics pattern that inflates very quickly and basically self-similarly after the first ∼100 ms.

4.2.3. Ejecta Anisotropy and Kick Systematics

Aside from listing the absolute values of the final NS velocity, v_ns, and acceleration, a_ns, at the end of our simulations, Table 3 also compiles quantities that are diagnostically relevant for the SN explosion, i.e., the explosion energy, E_exp; the time for the onset of the explosion, t_exp, which is defined as the instant when the (diagnostic) explosion energy exceeds 10⁴⁹ erg (the reader is reminded here that the corresponding post-bounce time is t_exp+18 ms); the baryonic mass of the NS, M_ns; the ejecta mass, M_ej, defined as the gas mass enclosed in the shell between the NS on the one side and the SN shock on the other. All of these quantities (except for t_exp) are given at the end of the simulations, which is after 1 s of computed evolution except for model O16-3D, which terminated after 750 ms. The values of v_ns result from the requirement of momentum conservation of ejecta gas plus NS in the hydrodynamical simulation (Equation (2)); the corresponding time derivative yields a_ns. Note that the notation $\langle X\rangle$ indicates that the tabulated value of quantity X is the average over the last 20 ms of the respective simulation.

The last two columns of Table 3 provide values for the quantities P_ej and α_ej, which are the (radial) ejecta momentum and the momentum-asymmetry parameter of the ejecta, respectively. These two quantities allow for a better understanding of the origin of systematic trends in the obtained NS kicks. Following Scheck et al. (2006), they are defined as

$$\begin{eqnarray}&&{\alpha }_{\mathrm{ej}}=\displaystyle \frac{| {{\boldsymbol{P}}}_{\mathrm{gas}}| }{{P}_{\mathrm{ej}}},\end{eqnarray} \tag{ 5 }$$

which relates the net gas momentum $| {{\boldsymbol{P}}}_{\mathrm{gas}}|$ (see Equation (2) and text below this equation) to the scalar quantity

$$\begin{eqnarray}&&{P}_{\mathrm{ej}}={\int }_{{R}_{\mathrm{ns}}}^{{R}_{{\rm{s}}}}\rho | {\boldsymbol{v}}| \,{dV}.\end{eqnarray} \tag{ 6 }$$

According to Equation (6), P_ej represents the total (basically radial, because $| {\boldsymbol{v}}| ={v}_{r}$ to high accuracy) momentum available in the volume between NS surface (at R_ns(θ, ϕ)) and SN shock (at radius R_s(θ, ϕ)). Using these definitions, we can write the magnitude of the NS kick velocity in terms of the radial ejecta momentum and the momentum-asymmetry parameter as

$$\begin{eqnarray}&&| {{\boldsymbol{v}}}_{\mathrm{ns}}| ={\alpha }_{\mathrm{ej}}\,\displaystyle \frac{{P}_{\mathrm{ej}}}{{M}_{\mathrm{ns}}}.\end{eqnarray} \tag{ 7 }$$

The parameter α_ej (which can vary between 0 and 1) is suitable to characterize the ejecta morphology. Under the assumption that the explosion energy is dominated by the kinetic energy of the expanding (anisotropic) ejecta¹¹ (or just a constant fraction of it), i.e., ${E}_{\exp }\propto {E}_{\mathrm{kin}}$, P_ej is related to the explosion energy by ${E}_{\exp }\propto \tfrac{1}{2}{P}_{\mathrm{ej}}^{2}{M}_{\mathrm{ej}}^{-1}$. This implies

$$\begin{eqnarray}&&{P}_{\mathrm{ej}}\propto \sqrt{2\,{M}_{\mathrm{ej}}\,{E}_{\exp }}.\end{eqnarray} \tag{ 8 }$$

The excellent validity of this relation for ECSNe can be concluded from the discussion following below and the results shown in Figure 9. For more massive Fe-core SN explosions, it will be demonstrated in a forthcoming paper (A. Wongwathanarat & H.-T. Janka 2018, in preparation).

Zoom In Zoom Out Reset image size

In neutrino-driven explosions, neutrino heating lifts matter from the vicinity of the PNS to a gravitationally less or even marginally bound state (i.e., to a specific binding energy closer to zero), while subsequent energy release by the recombination of free nucleons to heavy nuclei (up to ∼8.8 MeV per nucleon) provides the excess energy for the explosion (see Janka 2001; Marek & Janka 2009; Müller 2015).¹² Therefore, one would expect that the explosion energy grows with the mass of the neutrino-heated ejecta,

$$\begin{eqnarray}&&{E}_{\exp }\propto {M}_{\mathrm{ej}}\sim {M}_{{\rm{s}}}-{M}_{\mathrm{ns}},\end{eqnarray} \tag{ 9 }$$

where M_s is the mass enclosed by the SN shock during the first second(s) of the explosion. For the considered ECSN simulations, this expectation is confirmed by Figure 9 (lower left panel), which shows that E_exp increases linearly with decreasing NS mass M_ns. As more matter is expelled in the explosion by neutrino heating instead of being integrated into the NS, the explosion energy increases by the same degree. The extra mass expelled enlarges the mass available between NS and shock, M_s − M_ns.

Using the proportionality relation of Equation (9) in (8) yields

$$\begin{eqnarray}&&{P}_{\mathrm{ej}}\propto {E}_{\exp }.\end{eqnarray} \tag{ 10 }$$

This relation is fully in line with Figure 9, lower right panel. Equation (10) in (7) finally leads to

$$\begin{eqnarray}&&| {{\boldsymbol{v}}}_{\mathrm{ns}}| \propto \displaystyle \frac{{\alpha }_{\mathrm{ej}}}{{M}_{\mathrm{ns}}}\,{E}_{\exp };\end{eqnarray} \tag{ 11 }$$

see also Equation (11) of Janka (2017b).

Our result for Equation (11) explains what we see for our large sample of 2D models in the upper left panel of Figure 9: Stochastic variations of $| {{\boldsymbol{v}}}_{\mathrm{ns}}|$ stem from such variations in the asymmetry parameter α_ej (Figure 9, upper right panel), on top of which a linear trend of increase of $| {{\boldsymbol{v}}}_{\mathrm{ns}}|$ with E_exp is superimposed. Because of the tendency toward higher kicks at higher explosion energies, the variance of the kicks also exhibits this trend.

For all five energy sets, the kick velocity of the 3D model is below the average of the 2D models. This can be expected from the axisymmetric geometry of 2D simulations, which favors the formation of downflows at the poles (see Figure 2). In that sense, O8-3D is peculiar. Because of the stochastic development of ejecta asymmetries in the nonlinear stage of the growth of hydrodynamic instabilities, model O8-3D accidentally is the only 3D case that features a downflow at the end of the simulation that is similarly prominent to those in many of the 2D models. For this reason, model O8-3D is the only 3D case with an NS kick close to the mean value of the 2D runs.

The values of the asymmetry parameter α_ej (Figure 9, upper right panel) are below ∼3.5%; on average, they are around 1%. This is a factor of 10 smaller than those in typical iron-core SN explosions, where values up to roughly ∼35% were found, with an average around 10%; see Figure 11 of Scheck et al. (2006).

An interesting question concerns a possible dependence of the values of the momentum-asymmetry parameter α_ej on the explosion energy. The asymmetries of the 2D models exhibit a (slight) indication of a maximum at intermediate explosion energies. The 3D models seem to mirror such a trend, but this has to be taken with caution in view of the fact that we have only one 3D simulation for each value of the explosion energy. Moderate E_exp therefore might favor the highest anisotropies. Such a local peak could be imagined as a consequence of two competing effects: First, weak neutrino heating for cases with low explosion energies hampers the growth of large asymmetries by strong convective activity. In contrast, high explosion energies due to strong neutrino-energy deposition imply very fast expansion of the convectively perturbed, neutrino-heated ejecta. Therefore, the matter in the close vicinity of the NS is quickly replaced by the spherically symmetric neutrino-driven wind, which does not contribute to the recoil acceleration of the NS. Hence, O3 models (i.e., models with the lowest explosion energies) yield NS kicks $| {{\boldsymbol{v}}}_{\mathrm{ns}}|$ that are small, on average, because of low values of α_ej with little variance as well as low available momentum P_ej. Likewise, O12 and O16 models (with their highest explosion energies in our sample) feature relatively small α_ej, but their average NS kicks are on the larger side, due to the higher values of P_ej. Furthermore, the variance of the kick velocities is considerable (upper left panel of Figure 9). This means that increasing ejecta momentum compensates for the decrease of the anisotropy parameter to produce relatively large kicks. Of course, conclusive results for the dependence of α_ej on E_exp would require still more models—in particular, an extensive sample of 3D calculations.

What caused the kick of the Crab pulsar? In view of our results presented in Section 4.2, the answer is far from being clear. If Crab is the remnant of an O–Ne–Mg core progenitor, a hydrodynamical kick mechanism for its pulsar is strongly disfavored by the extremely low NS recoil velocities obtained in our simulations.

A possible alternative might be a postnatal acceleration of the pulsar due to electromagnetic radiation produced by an off-centered, rotating magnetic dipole (Harrison & Tademaru 1975). The NS kick by this electromagnetic rocket effect is attained on the spin-down timescale and is imparted along the spin axis of the pulsar. The NS obtains its proper motion at the expense of kinetic energy associated with its rotation. This scenario appears attractive because of the impressive radiation produced by the Crab pulsar even at the present age, and because of the observed rough spin-kick alignment of this compact object, provided this alignment is not purely incidental (for a discussion with references, see Lai et al. (2001)).

However, a quantitative analysis reveals that this explanation is unlikely. Pushing all involved parameters to their extreme limits, a pulsar kick of at most

$$\begin{eqnarray}&&{v}_{\mathrm{HT},\mathrm{kick}}^{\max }\approx 5.05\,{R}_{12}^{2}{\left(\displaystyle \frac{20\mathrm{ms}}{{T}_{\mathrm{ns},{\rm{i}}}}\right)}^{2}\,\displaystyle \frac{\mathrm{km}}{{\rm{s}}}\end{eqnarray} \tag{ 12 }$$

can be expected by the Harrison–Tademaru effect (Lai et al. (2001), Equation (5), with pulsar spin frequency $\nu ={T}_{\mathrm{ns}}^{-1}$ and $\nu \ll {\nu }_{{\rm{i}}}$). In Equation (12), R₁₂ = R_ns/(12 km) is the NS radius and T_ns,i is the initial NS spin period. This means that, in order to achieve a recoil of 160 km s⁻¹ via the emission of electromagnetic energy, an initial value of T_ns,i ∼ 3.5 ms is required, in conflict with the estimated ∼19 ms birth period of the Crab pulsar (Manchester & Taylor 1977; Bejger & Haensel 2003). Moreover, the output of spin-down energy of a 3.5 ms pulsar, ${\rm{\Delta }}{E}_{\mathrm{rot}}\,\sim$ $\displaystyle \frac{1}{2}{I}_{\mathrm{ns}}(2\pi$/${T}_{\mathrm{ns},{\rm{i}}}{)}^{2}\approx \tfrac{1}{2}$ $\left(\tfrac{2}{5}{M}_{\mathrm{ns}}{R}_{\mathrm{ns}}^{2}\right)$ ${(2\pi /{T}_{\mathrm{ns},{\rm{i}}})}^{2}$ $\sim 2.8\times {10}^{51}{M}_{1.5}{R}_{12}^{2}$ ${({T}_{\mathrm{ns},{\rm{i}}}/3.5\mathrm{ms})}^{-2}$ erg (with the NS mass M_1.5 = M_ns/(1.5 M_⊙) and NS moment of inertia I_ns), is incompatible with the energetics and dynamics of the Crab Nebula (Yang & Chevalier (2015); but see Atoyan (1999) for a different opinion, arguing that basically all of the injected energy could have been radiated away).

A second alternative possibility for explaining the proper motion of the Crab pulsar could be the SN-induced breakup of a close binary system, in which the Crab progenitor would have been one component. This Blaauw effect (Blaauw 1961) could account for the inferred ∼160 km s⁻¹ of the NS, but evidence for a former companion star does not exist. Moreover, if this scenario were true, the current interpretation of the observed southeast-northwest asymmetry of the Crab's synchrotron nebula as a consequence of the proper motion of the pulsar would require revision. A displacement of the pulsar from the explosion center (Hester 2008, Figure 7 therein) would not be compatible with a binary disruption by the Crab SN, in which case the NS and the center of the expanding ejecta cloud would move with the same speed as the progenitor's previous orbital motion in the binary system.

A third possible alternative could be a low-mass iron-core progenitor for the origin of Crab. SN explosions of such progenitors in the ∼9 M_⊙ to ∼12 M_⊙ range (Woosley & Heger 2015) eject similar amounts of stellar debris with similarly low explosion energies (of order ∼10⁵⁰ erg; see Melson et al. (2015b), Radice et al. (2017)) as ECSNe, and they also underproduce ⁵⁶Ni compared to core-collapse SNe of more massive progenitors (Ertl et al. 2016; Sukhbold et al. 2016). Nevertheless, low-mass iron-core progenitors can possess considerably higher values for the core compactness (i.e., shallower density profiles exterior to the iron core). There is a considerable variation in the density profiles between different representatives in this mass range (Sukhbold et al. 2016, Figure 3; Janka 2017a, Figure 2), corresponding to an increase of ξ_1.5 from about 7.5 × 10⁻³ for a 9.0 M_⊙ progenitor to ξ_1.5 ∼ 0.094 for a 9.5 M_⊙ star and ξ_1.5 between 0.44 and 0.58 for models between 10.0 M_⊙ and 12.0 M_⊙. Therefore, such stars might explode with sufficient ejecta momentum asymmetry to produce NS kick velocities in the ballpark of 160 km s⁻¹. This argument receives support from the results of Suwa et al. (2015) and Müller et al. (2018), who reported (non-converged) kick velocities of several 10 km s⁻¹ up to ∼75 km s⁻¹ within less than one second of post-bounce evolution for a small set of 2D neutrino-driven SN simulations of ultra-stripped SNe Ic progenitors (with compactness values ξ_1.4 between ∼0.47 and ∼0.97). However, further investigation is needed to clarify whether low-mass iron-core progenitors are compatible with the ejecta composition of the Crab Nebula, whose high helium abundance and relatively oxygen-poor filaments served as an argument for an O–Ne–Mg-core progenitor in Nomoto et al. (1982) and Davidson & Fesen ( 1985); see also the concise overview of Crab properties and corresponding references provided by Smith (2013).

Finally, a fourth possibility is that the Crab pulsar attained its kick via anisotropic neutrino emission during the neutrino-cooling phase of the hot PNS. The corresponding neutrino-induced kick velocity can be written as

$$\begin{eqnarray}&&{v}_{\nu ,\mathrm{kick}}=167\,\displaystyle \frac{{\alpha }_{\nu }}{0.005}\,\displaystyle \frac{{E}_{\nu }}{3\times {10}^{53}\,\mathrm{erg}}{\left(\displaystyle \frac{{M}_{\mathrm{ns}}}{1.5{M}_{\odot }}\right)}^{-1}\,\displaystyle \frac{\mathrm{km}}{{\rm{s}}},\end{eqnarray} \tag{ 13 }$$

where α_ν = P_z,ν/P_ν is the neutrino momentum asymmetry parameter and P_ν = E_ν/c is the total radial momentum associated with a neutrino-energy loss E_ν (Scheck et al. 2006). The parameter α_ν is related to the amplitude a_d of the dipole component (oriented in the z-direction) of a neutrino emission asymmetry, ${{dL}}_{\nu }(t)/d{\rm{\Omega }}={L}_{\nu }(t)(1+{a}_{{\rm{d}}}\cos \theta )/(4\pi )$, by the relation α_ν = a_d/3. Therefore, Equation (13) implies that an emission dipole of roughly 1.5% of the monopole, applied to the total neutrino-energy loss E_ν, is sufficient to obtain a kick of the magnitude of the Crab pulsar.

One way to generate considerable dipolar asymmetry of the neutrino emission is by the presence of very strong magnetic fields in the NS interior (e.g., Bisnovatyi-Kogan 1993; Socrates et al. 2005; Kusenko et al. 2008; Sagert & Schaffner-Bielich 2008; Maruyama et al. 2012), which can modify the neutrino interactions in the dense matter. Corresponding neutrino-induced kicks of some 100 km s⁻¹ typically require dipolar field asymmetries of order 10¹⁵ G or more. However, this requirement is in conflict with the fact that the Crab pulsar has an ordinary dipolar surface field of about 4 × 10¹² G (e.g., Michel 1991), and interior fields of order 10¹⁵ G would be purely speculative. On the other hand, neutrino-emission dipoles of several percent amplitude have been found in association with the newly discovered dipolar lepton-emission self-sustained asymmetry (LESA); see Tamborra et al. (2014a) and Janka et al. (2016), as well as the recent confirmation by O'Connor & Couch (2018) in full-scale 3D simulations of NS birth. LESA does not require the presence of a magnetic field. The corresponding emission dipole exhibits a stable or only slowly changing direction and can amount to a few percent of the added luminosities of neutrinos and antineutrinos of all species. However, current 3D simulations have been able to follow this phenomenon in the nascent NS only for the first 500–700 ms after core bounce. It is therefore unclear how long a significant dipolar asymmetry of the neutrino emission can persist and whether it can reach an overall amplitude on the order of 1%–2% of the total (i.e., time-integrated) neutrino-energy loss, which is needed to explain an NS kick of about 160 km s⁻¹ (see Equation (13)). If LESA is a generic phenomenon during the neutrino-cooling phase of all hot, newborn NSs, it would lead to a relatively high "floor value" for NS kicks; see Bray & Eldridge (2016) and Bray & Eldridge (2018) for a corresponding conjecture based on combining measured NS proper motions and stellar population modeling of stars evolving to SNe. However, such a ubiquitous NS acceleration mechanism, which would act in all cases, might be in conflict with constraints associated with a population of low-kick NSs in binary systems (Podsiadlowski et al. 2004; Schwab et al. 2010; Beniamini & Piran 2016; Tauris et al. 2017; Kruckow et al. 2018).

In summary, the origin of the kick of the Crab pulsar remains a puzzle and could be a smoking gun for either the nature of the progenitor of the Crab SN or neutrino-related emission asymmetries of the nascent NS.

We have presented results of 40 2D and five 3D simulations of neutrino-driven ECSN explosions of an O–Ne–Mg-core progenitor, considering models with parametrically tuned explosion energies between about 3 × 10⁴⁹ erg and roughly 1.6 × 10⁵⁰ erg. This range of energies brackets the results of previous self-consistent core-collapse and explosion simulations (Kitaura et al. 2006; Janka et al. 2008; Fischer et al. 2010; Hüdepohl et al. 2010; von Groote 2014; Radice et al. 2017) and the estimated explosion energy of the Crab SN (Yang & Chevalier 2015), which is discussed as a possible candidate for an ECSN (e.g., Hillebrandt 1982; Nomoto et al. 1982; Smith 2013; Tominaga et al. 2013).

The goal of our study was to explore the hydrodynamic NS kicks that are associated with anisotropic mass ejection in this type of SN explosion. The asymmetries originate from a brief episode of growth of Rayleigh–Taylor instability in the convectively unstable neutrino-heated region exterior to the gain radius. The emerging inhomogeneities in the mass and energy distributions of the neutrino-heated ejecta exert hydrodynamical and gravitational forces on the NS and thus accelerate the NS via the gravitational tugboat mechanism. The typical NS kick velocities that we obtained in our ECSN models are on the order of kilometers per second.

Our results support previous theoretical arguments that NS kicks in ECSNe should be smaller on average, compared to NS kicks associated with SNe of more massive iron-core progenitors (Podsiadlowski et al. 2004; Janka 2017b). The extremely low compactness values of the stellar core (ξ_1.4 = 1.1 × 10⁻⁵ and ξ_1.5 = 6.6 × 10⁻⁶) enable a very rapid outward acceleration of the SN shock, which also allows the postshock layer to expand rapidly in the wake of the shock. For this reason, the merging process of the Rayleigh–Taylor plumes of neutrino-heated, buoyant plasma to larger structures freezes out in a higher-order spherical harmonics pattern with insignificant amplitudes of dipolar or quadrupolar deformation modes.

This ECSN-typical dynamical scenario has three consequences that conspire to disfavor large NS kicks: First, the dipolar momentum asymmetry (measured by the parameter α_ej; Equation (5)) is only on the order of one percent, and thus about 10 times smaller than for typical iron-core SNe. Second, the low explosion energy of ECSNe (around 10⁵⁰ erg) implies a low radial momentum of the ejecta and also reduces the NS kick (see Equations (10) and (11)). Third, the extremely fast expansion of the convectively perturbed ejecta diminishes the time over which the anisotropic forces responsible for the gravitational tugboat acceleration can act on the NS. These effects explain why our simulations yield NS kick velocities of only a few kilometers per second at most, which is a factor of 100 below those typical of iron-core SNe (Scheck et al. 2006; Wongwathanarat et al. 2013). Our extensive set of 40 2D and five 3D simulations for a relevant range of explosion energies accounts for stochastic variations of the ejecta asymmetry, and thus renders this finding a statistically consolidated result. In view of our pool of results, testing different explosion energies and initial perturbations, it appears very unlikely that current ECSN models can explain NS kicks of more than 100 km s⁻¹ with any significant statistical probability.

The low-velocity NS kicks that we obtained in our ECSN models are compatible with the possibility of a dichotomous distribution of NS kicks (e.g., Katz 1975, 1983; Podsiadlowski et al. 2004; Schwab et al. 2010; Beniamini & Piran 2016, and references therein), and they also back up the possibility of a large population of double compact objects as gravitational-wave sources (e.g., Tauris et al. 2017; Kruckow et al. 2018). However, such low NS kicks are in tension with the hypothetical explanation of the Crab Nebula as the gaseous remnant of an ECSN (e.g., Hillebrandt 1982; Nomoto et al. 1982; Smith 2013; Tominaga et al. 2013), because the Crab pulsar's measured proper motion of ∼160 km s⁻¹ is two orders of magnitude higher than the average NS kicks found in our ECSN explosion models. Such velocities are out of reach by hydrodynamical kicks in ECSNe even for statistical outliers.

A postnatal acceleration of the Crab pulsar by the electromagnetic rocket effect (Harrison & Tademaru 1975) might appear attractive because of the apparent alignment of the pulsar's spin and kick. However, this explanation is strongly disfavored by its requirement of a short initial spin period of less than 4 ms (see detailed discussion in Section 5), whose associated release of electromagnetic power is more than an order of magnitude higher than limits set by a detailed energetic analysis of the Crab Nebula (Hester 2008; Yang & Chevalier 2015).

Alternatively, attributing the proper motion of the Crab pulsar to the disruption of a binary system in the event of SN 1054 is in tension with the current picture of the internal structure and dynamics of the Crab Nebula. If the NS's space velocity had originated from the disruption of a close binary system, the pulsar would not move relative to the mass-center of the ejecta gas, and a displacement of the pulsar out of the center of the explosion (Hester 2008, Figure 7) would not apply. Moreover, the spin-kick alignment of the Crab pulsar would imply that the spin axis was close to the orbital plane, which is not the most plausible configuration of spin and orbital angular momentum.

Two other possible explanations of the Crab pulsar's proper motion appear more likely. Either the progenitor of SN 1054 was a low-mass star with an iron core instead of the hypothesized O–Ne–Mg-core progenitor, or the Crab pulsar received its kick velocity via anisotropic neutrino emission rather than being kicked by the hydrodynamical mechanism associated with asymmetric mass ejection in the SN explosion.

Neutrino-induced NS kicks of ≳100 km s⁻¹ may be a consequence of anisotropic neutrino emission associated with an ultra-strong dipole magnetic field of about 10¹⁵ G or more (e.g., Bisnovatyi-Kogan 1993; Socrates et al. 2005; Kusenko et al. 2008; Sagert & Schaffner-Bielich 2008; Maruyama et al. 2012). Because the Crab pulsar is not a magnetar, the existence of such enormous magnetic fields in its interior would be pure speculation. A dipolar asymmetry of the neutrino emission is also associated with the LESA instability newly discovered in 3D SN simulations (Tamborra et al. 2014a, 2014b; Janka et al. 2016; O'Connor & Couch 2018), which does not depend on the presence of a strong magnetic field. If this asymmetry persists for seconds (and thus much longer than the ∼0.5 s over which it was feasible to track the PNS evolution in recent 3D models), it is conceivable that LESA could well account for NS kick velocities of 100–200 km s⁻¹. It may be speculated whether this effect could be connected to a constant floor value of the NS kick velocity of more than 100 km s⁻¹ that has been surmised by Bray & Eldridge (2016). An open question remains as to whether such a substantial floor value by a neutrino-induced kick would be compatible with the large population of NSs in globular clusters (Katz 1975, 1983; Pfahl et al. 2002a; Kuranov & Postnov 2006) and with the dichotomy of kick velocities advocated by Katz (1975), Pfahl et al. (2002b), Podsiadlowski et al. (2004), Beniamini & Piran (2016), and Tauris et al. (2017), who invoke a low-kick population of NSs to explain the existence of a certain class of high-mass X-ray binaries and the properties of binary NS systems.

Solar-metallicity iron-core progenitors between ∼9 M_⊙ and ∼12 M_⊙ possess much shallower density profiles exterior to their iron cores, and therefore have higher core-compactness values than O–Ne–Mg-core progenitors. Ergo, it is conceivable that some of these cases might explode with sufficiently large anisotropy of the mass ejection to account for NS kicks beyond 100 km s⁻¹. Of course, referring to the possibility of sufficiently large NS kick velocities in low-mass iron-core SNe requires confirmation by an extensive set of 3D simulations of such explosions. Moreover, a better understanding is needed whether the explosion properties of SN 1054 (energy, iron production, light curve) and the chemical composition of the Crab Nebula are compatible with an origin from a low-mass iron-core progenitor.

In summary, our work reveals problems with the widely-used interpretation of the Crab Nebula as a remnant of an ECSN. The Crab pulsar's proper motion cannot be explained by hydrodynamic kicks associated with the ejecta asymmetries of ECSNe. Considering the pulsar velocity as a consequence of a natal kick, a connection of Crab to the explosion of a low-mass iron-core progenitor seems more plausible. Such SNe can also have low explosion energies and correspondingly little nickel production (Janka et al. 2012; Melson et al. 2015b; Müller 2016; Radice et al. 2017; Wanajo et al. 2018), compatible with recent reexamination of the evolution of the Crab pulsar wind nebula and its interaction with the SN ejecta (Yang & Chevalier 2015), which corroborates the notion that SN 1054 was a low-energy explosion (∼10⁵⁰ erg), and that significant energy in extended fast material around the Crab is unlikely. To assess the conflicting progenitor possibilities, further studies of the Crab Nebula by detailed observations, particularly regarding the chemical composition, are desirable. Moreover, self-consistent 3D SN models are needed for a large sample of low-mass (O–Ne–Mg and Fe-core) progenitors, clarifying the question of how steep the increase of the maximum achievable NS kick velocity with higher core-compactness can be.

We are grateful to Annop Wongathanarat for providing the 3D Yin-Yang version of the Prometheus-HOTB code, to Ken Nomoto for providing the stellar progenitor model, to Lorenz Hüdepohl for the post-bounce data of the collapsed progenitor star, to Thomas Ertl and Tobias Melson for helpful discussions and valuable information, and to Paz Beniamini, Jonathan Katz, Bernhard Müller, and Alexander Tutukov for insightful comments on the first arXiv version. We further thank Elena Erastova (Max Planck Computing and Data Facility) for the 3D visualization of the simulation shown in Figure 5. This research was supported by the Deutsche Forschungsgemeinschaft through Excellence Cluster Universe (EXC 153; http://www.universe-cluster.de/) and Sonderforschungsbereich SFB 1258 "Neutrinos and Dark Matter in Astro- and Particle Physics," and by the European Research Council through grant ERC-AdG No. 341157-COCO2CASA. The computations of the presented models were carried out on Hydra of the Max Planck Computing and Data Facility (MPCDF) Garching.

Software: Prometheus-HOTB (Fryxell et al. 1989; Janka & Müller 1996; Kifonidis et al. 2003; Scheck et al. 2006; Arcones et al. 2007; Wongwathanarat et al. 2013), VisIt (https://wci.llnl.gov/simulation/computer-codes/visit/).