Three-dimensional Models of Core-collapse Supernovae From Low-mass Progenitors With Implications for Crab

We present 3D full-sphere simulations of supernovae of non-rotating low-mass (~9 Msun) progenitors, covering the entire evolution from core collapse through bounce and shock revival, through shock breakout from the stellar surface, until fallback is completed several days later. We obtain low-energy explosions [~(0.5-1.0)x 10^{50} erg] of iron-core progenitors at the low-mass end of the core-collapse supernova (LMCCSN) domain and compare to a super-AGB (sAGB) progenitor with an oxygen-neon-magnesium core that collapses and explodes as electron-capture supernova (ECSN). The onset of the explosion in the LMCCSN models is modelled self-consistently using the Vertex-Prometheus code, whereas the ECSN explosion is modelled using parametric neutrino transport in the Prometheus-HOTB code, choosing different explosion energies in the range of previous self-consistent models. The sAGB and LMCCSN progenitors that share structural similarities have almost spherical explosions with little metal mixing into the hydrogen envelope. A LMCCSN with less 2nd dredge-up results in a highly asymmetric explosion. It shows efficient mixing and dramatic shock deceleration in the extended hydrogen envelope. Both properties allow fast nickel plumes to catch up with the shock, leading to extreme shock deformation and aspherical shock breakout. Fallback masses of<~5x10^{-3} Msun have no significant effects on the neutron star (NS) masses and kicks. The anisotropic fallback carries considerable angular momentum, however, and determines the spin of the newly-born NS. The LMCCSNe model with less 2nd dredge-up results in a hydrodynamic and neutrino-induced NS kick of>40 km/s and a NS spin period of ~30 ms, both not largely different from those of the Crab pulsar at birth.

The explosion is powered by gravitational energy, which is released when the core of the star collapses to a compact remnant (a neutron star or black hole), and a fraction of which is transferred to the ejecta by neutrino-energy deposition (Bethe & Wilson 1985;Colgate & White 1966). In the past six decades, numerous studies have focused on the collapse phase and subsequent evolution using 1D simulations. Over the past three decades, multi-dimensional simulations with successively improved treatment of the microphysics have driven our understanding of the explosion mechanism (see e.g. Janka et al. 2012;Janka 2012;Burrows 2013;Janka et al. 2016;Müller 2016). A close analysis of these models led to the discovery of new hydrodynamic instabilities such as the standing accretion shock instability (SASI) (Blondin et al. 2003;Blondin & Mezzacappa 2007;Ohnishi et al. 2006;Foglizzo et al. 2007;Fernández 2010). The increase of computational capabilities in the recent years along with new developments for neutrino transport methods (see, e.g. Buras et al. 2006;Takiwaki et al. 2012;O'Connor & Couch 2018a;Skinner et al. 2019;Glas et al. 2019a and references therein) have enabled full 3D simulations of the early explosion phase (see e.g. Takiwaki et al. 2014;Couch & O'Connor 2014;Melson et al. 2015a,b;Lentz et al. 2015;Summa et al. 2016;Roberts et al. 2016;Müller et al. 2017Müller et al. , 2018Vartanyan et al. 2018;Ott et al. 2018;O'Connor & Couch 2018b;Melson et al. 2019;Vartanyan et al. 2019;Burrows et al. 2019).
Motivated by the historical SN1987A and its progenitor detection, a variety of studies in two and three dimensions also investigated the propagation of the shock wave from its initiation to its breakout from the stellar surface in 15-20 M blue supergiant (BSG) models, which are suitable as progenitors of SN1987A. In first studies, the blast wave was launched through artificial energy deposition near the center (see e.g. Nomoto et al. 1987Nomoto et al. , 1988Arnett et al. 1989b;Müller et al. 1991b;Fryxell et al. 1991), later the explosion was initiated with the neutrino-driven mechanism (see Kifonidis et al. 2003;Hammer et al. 2010;Wongwathanarat et al. 2013;Wongwathanarat et al. 2015). Müller et al. (2018) followed the long-time evolution of the explosion of ultra-stripped progenitors by 3D simulations, motivated by the importance of such stars in understanding the progenitor systems of the recent detections of NS-NS mergers by LIGO/Virgo (GW170817, Abbott et al. 2017, andGW190425, Ligo Scientific Collaboration &VIRGO Collaboration 2019).
These theoretical works showed that supernova explosions are by far not spherical events as previously thought. Three-dimensional instabilities facilitate the explosion (Herant et al. 1994;Burrows et al. 1995;Janka & Müller 1996) and are a necessary ingredient to explain the clumpiness and mixing found in photospheric emission (Utrobin et al. 2015(Utrobin et al. , 2017 and spectral analyses of the nebular phase of core-collapse events (Jerkstrand 2017). Wongwathanarat et al. (2015) showed that the final ejecta distribution carries imprints of the asphericities produced during the onset of the explosion (t ∼ 1 s), which are further modified during later phases. Depending on the detailed progenitor structure, hydrodynamic instabilities arising at the composition interfaces, such as the Rayleigh-Taylor (RT) instability and the Richtmeyer-Meshkov (RM) instability, shape the final spatial and velocity distributions of nucleosynthetic products. The resulting ejecta morphology ranges from quasi-spherical ejecta  to strongly pronounced RT-fingers including cases that resemble the geometry found in Cas A Grefenstette et al. 2017). Due to the highly nonlinear and stochastic behaviour of non-radial instabilities and turbulence during the onset of the explosion phase and the subsequent evolution of RT instabilities, which depend on the progenitor structure, a clear connection between the asymmetries, and thus the degree of mixing, and the progenitor properties has still to be worked out.
In this paper, we consider CCSN progenitors with initial masses near the low-mass end of about 9-10 M , where around 20% of all CCSNe are thought to occur (assuming a Salpeter initial mass function Salpeter 1955 and an upper mass limit of ∼20 M for CCSNe), to study the differences in their development of mixing instabilities during the explosion. To this end, we compare ECSNe from a super-AGB progenitor with an ONeMg core and CCSNe from red supergiants (RSG) progenitors with iron cores, all in a zero-age main sequence mass range around 9 M .
The evolution of stars with masses 12 M is very sensitive to the initial stellar mass, various pulsational instabilities, and massloss phenomena (Woosley & Heger 2015). Iron-core CCSN progenitors with initial masses around 9-10 M ignite oxygen burning off-center in contrast to their more massive (M > 15 M ) counterparts. After oxygen burning, silicon ignites in a degenerate flash which might, in some cases, lead to additional mass loss in the last decade of evolution or is speculated to even eject parts of the hydrogen envelope (Woosley & Heger 2015).
The 8.8 M progenitor of Nomoto (1984) is even more peculiar. It experiences several thermal pulses and off-center ignition of fusion material. In the end, it has a degenerate ONeMg-core surrounded by a dilute and extended hydrogen envelope.
When the core approaches its Chandrasekhar mass, electron captures on 24 Mg and 20 Ne via the reaction chains 24 Mg(e − , ν e ) 24 Na(e − , ν e ) 24 Ne and 20 Ne(e − , ν e ) 20 F(e − , ν e ) 20 O (Miyaji et al. 1980) destabilize the core due to a reduction of the effective adiabatic index of the electron-degeneracy dominated gas pressure. Continuous electron capture on 20 Ne, which further reduces the pressure support, works against the now beginning oxygen-burning as temperatures increase during the collapse. Simulations in 1D (see Kitaura et al. 2006;Hüdepohl et al. 2010;Fischer et al. 2010) and 2D (see Janka et al. 2008;Radice et al. 2017) suggest that the collapse proceeds despite the oxygen burning. Jones et al. (2016) simulated the deflagration of oxygen in ONeMg cores with different core densities. At log 10 (ρ c /g cm −3 ) = 9.95 and lower densities their cores do not collapse but get partly unbound due to the inefficient semi-convective mixing during the electron-capture phase and the resulting strong thermonuclear runaway. Only when the central densities are higher than this threshold value of ρ c , the core is found to collapse to a proto-neutron star (PNS). Recently, Kirsebom et al. (2019) investigated the influence of a newly measured strong transition between the ground states of 20 Ne and 20 F on the electron-capture rate and thus on the evolution of ONeMg cores. Adding the new transition increases the likelihood that the star is (partially) disrupted by a thermonuclear explosion (termed tECSN) rather than collapsing to form a PNS. However, Zha et al. (2019), using state-of-the-art electron-capture rates including the latest rate for the second forbidden transition of 20 Ne(e − , ν e ) 20 F from Suzuki et al. (2019), found that the oxygen deflagration starting from log 10 (ρ c /g cm −3 ) > 10.01 (< 10.01) leads to collapse (thermonuclear explosion). Their estimate of the central density when the oxygen deflagration is initiated in an evolving ONeMg core exceeds this critical value. Therefore they conclude that ONeMg cores are likely to collapse.
For this reason, in our study we assume that the ONeMg core collapses to a PNS, leading to a "collapse ECSN" (cECSN). This assumption receives additional motivation by the fact that recent studies considering the galactic chemical evolution of the Milky Way stress the importance of cECSNe to reproduce the solar abundances of several important and problematic isotopes including, e.g., 48 Ca, 50 Ti, and several of the isotopes from Zn to Zr (Jones et al. 2019).
Despite the narrow range of central densities for which cEC-SNe are expected to occur (Leung et al. 2020), and despite the open questions associated with a variety of competing processes that decide about collapse or thermonuclear explosion and that depend strongly on many uncertain aspects of the employed physics, connections to cECSNe have been made for observations of SN 1994N, 1997D, 1999br, 1999eu, 2011dc and 2005cs (Stevenson 2014. However, comparisons of the nebular spectra of some of these cases with 1D neutrino-driven SN models are ambiguous or disfavor the link to cECSNe (Jerkstrand et al. 2018). Also SN1054 (the Crab) has been speculated to be a cECSN (Nomoto et al. 1982;Hillebrandt 1982;Tominaga et al. 2013;Smith 2013), although such an interpretation is in conflict with results by Gessner & Janka (2018) for the maximum kick velocity of PNSs produced by cEC-SNe.
With the help of full-sphere three-dimensional simulations we aim at investigating the following questions: • What are the differences in the early stages (first seconds) of the explosion in low-mass Fe-core and ONeMg-core progenitors?
• What is the influence of the different progenitor structures on the long-time evolution of the explosion? In particular, what is the influence on the formation of reverse shocks and the efficiency of outward mixing of neutrino-heated material?
• Are CCSNe of low-mass progenitors able to produce highly asymmetric ejecta and strong radial mixing of metals similar to findings for more massive RSG and BSG stars in previous studies?
• How do the properties of the compact remnants change on long time-scales due to the fallback of matter? Are there significant changes to the remnants' mass, kick, and angular momentum?
The structure and contents of our paper are the following: In Section 2, the basic properties of the considered models of non-rotating, low-mass (super-AGB and RSG) progenitors are introduced. Section 3 provides a brief description of the numerical methods and input physics used in our simulations. Section 4 contains our results for the first second(s) of the explosion, focusing on shock dynamics, explosion energies, neutrino emission, PNS properties and the chemical composition of the ejecta. In Section 5, for the first time in 3D explosion modeling, the SN evolution of low-mass super-AGB and RSG progenitors is described until and beyond shock breakout concerning the development of mixing instabilities and ejecta asymmetries, the spatial distribution of chemical elements, and the effects of fallback on the properties of the newly formed NSs. In Section 6, we briefly compare our results with previous studies, and in Section 7, we conclude with a summary and discussion. Several appendices contain basic information on more technical aspects concerning the simulation inputs and the analysis methods.

PROGENITORS
In this paper, we study an ECSN of a non-rotating 8.8 M super-AGB star (e8.8), which is constructed from the envelope model of Jones et al. (2013) and a collapsing core model (Leung et al. 2020;A. Tolstov, S.-C. Leung, and K. Nomoto, 2017, private communication), and two CCSNe resulting from non-rotating low-mass RSGs with iron cores (z9.6, s9.0), evolved to the onset of collapse by A. Heger (2015, private communication) and by Woosley & Heger (2015), respectively. The considered ECSN progenitor is explored here for the first time, whereas the explosions of the iron core progenitors were simulated in 3D with the V -P code and some results were published in Melson et al. (2015a) and Melson et al. (2019).

ONeMg-core progenitor
Model e8.8 is a solar-metallicity progenitor with a zero-agemain-sequence (ZAMS) mass of M ZAMS = 8.8 M . Its degenerate ONeMg-core undergoes electron capture which ignites the O-Ne deflagration at the center. The central density (ρ c ) and the electron fraction (Y e ), and thus the core mass M(ONeMg) at the ignition, depend on the convective stability criterion, extent of the convective mixing, and the electron capture rate . Leung et al. (2020) adopted the core with ρ c ∼ 10 9.975 g cm −3 , Y e = 0.496, and M(ONeMg) = 1.39 M , and calculated the propagation of the convective deflagration. The deflagration incinerates ONeMg-core material into nuclear statistical equilibrium (NSE). In the NSE region, electron capture on iron-group nuclei and free protons dominates nuclear energy release, thus inducing collapse. When the central density reaches ρ c ∼ 10 10.64 g cm −3 , the NSE region extends to ∼0.45 M . The SN simulations started from this progenitor condition. The structure of the core is well approximated by a spherical model.
When transferring the progenitor model of e8.8 to the SN simulations, the mass of the ONeMg-core was reduced to 1.34 M in the course of providing an easy-to-handle fit to the complex density structure. This has no relevant influence on the dynamical evolution of the explosion, as can be concluded from the close similarity of the explosion behavior of model e8.8 with that of a previous version of the progenitor (Nomoto 1984;Nomoto 1987;K. Nomoto, 2008, private communication), which is commonly termed e8.8 in the literature and e8.8 n in this publication. Model e8.8 n had a degenerate core of 1.375 M , but a smaller hydrogen envelope (see Figure 1). It is a reference case for simulations of ECSNe and was used in various studies, focusing on explosion properties (Kitaura et al. 2006;Fischer et al. 2010), nucleosynthesis (Janka et al. 2008;Wanajo et al. 2011;Wanajo et al. 2018), effects of microphysics (Hüdepohl et al. 2010;von Groote 2014;Radice et al. 2017), and PNS kicks (Gessner & Janka 2018). In order to compensate for possible uncertainties in our ECSN simulations with the new progenitor model of e8.8, we vary its explosion energy in a set of 2D simulations, denoted by e8.8 3 , e8.8 6 , e8.8 10 , and e8.8 15 for E exp = (3, 6, 10, 15) × 10 49 erg.
Progenitor model e8.8 has the very sharp density gradient at 1.34 M near the edge of its compact ONeMg-core (ξ 2.5 = 5.7 × 10 −6 , ξ 1.5 = 8.0 × 10 −6 ) 1 that is characteristic of such progenitors prior to the onset of core collapse. This steep gradient is a prominent feature of the density (ρ) profile in Figure 1. In the same same figure we also show the ρr 3 -profile as well as the electron fraction Y e and the temperature T. All of these quantities are displayed as functions of enclosed mass and radius. The inner ∼0.45 M of the degenerate core contain iron-group nuclei and α-particles in NSE.  Profiles of the temperaturei (T ), densityi (ρ), electron fraction (Y e ) and ρr 3 for the pre-collapse progenitor models as functions of radial coordinate (left panels) and mass coordinate (right panels). Indicated by dash-dotted, dashed, and dotted lines are the outer boundaries of the degenerate (iron or NSE), CO, and He cores, respectively. Note the huge differences in the density and ρr 3 -profiles between the progenitors with iron and ONeMg-cores, in particular just outside of the CO core. We show the difference in the core structures of our ONeMg-core models in a zoom of the ρ vs. M(r) profiles in the rightmost panel.
the NSE/ONeMg, CO/He, and He/H composition interfaces, respectively. We define the locations of the composition interfaces, similar to Wongwathanarat et al. (2015), as those positions at the bottom of the respective layers of the star where the mass fractions X i drop below half of their maximum values in the layer. The radial positions of the composition interfaces are summarized in Table 1. In the top panel of Figure 2 we present the composition of model e8.8, where we combine all elements with mass numbers greater than 28 into "iron-group" (IG) material or nuclei in NSE. The ONeMg-core is surrounded by very thin carbon and helium shells (M C ≈ 8.1×10 −3 M , M He ∼ 2.1×10 −6 M ) and a hydrogen (H+He) envelope (M H ≈ 4.49 M ). The total masses of the different nuclei present in the entire pre-SN model are listed in Table 5.
For collapse and post-bounce evolution of an ONeMg-core Pre-collapse composition of models e8.8 (top), z9.6 (middle) and s9.0 (bottom) as a function of radius (left) and enclosed mass (right). Note the broken horizontal axis of the right panels. We combine all chemical elements with mass numbers greater than 28 into the "iron-group" (IG). The dash-dotted, dashed, and dotted lines indicate the outer boundaries of degenerate (iron or NSE), CO, and He cores, respectively.
progenitor we employ in all of this paper the new progenitor model e8.8 in 1D, 2D, and 3D simulations with the P -HOTB code as detailed in Section 3.4. The profile of the old progenitor e8.8 n is shown in Figure 1 merely for illustration and reference.

Fe-core progenitors
As a second progenitor we employ a zero-metallicity M ZAMS = 9.6 M star, termed z9.6. It was first used by Janka et al. (2012) and was also considered in other studies such as Müller et al. (2013); Radice et al. (2017); Müller et al. (2019). This iron core progenitor is structurally similar to the ECSN model. It also shows a sharp decline of the density at the edge of its iron core, enabling low-energy explosions in 1D (Melson et al. 2015a;Radice et al. 2017). Evolved by A. Heger (2012, private communication) as an extension to Heger & Woosley (2010), the pre-SN model develops an iron core of about 1.30 M . The iron core is surrounded by a 0.061 M Si-layer, a 0.016 M CO-layer, and has a hydrogen-free helium layer of about 0.004 M below a 0.33 M convective H/Helayer with a hydrogen mass fraction of ∼6%, which is surrounded by a massive hydrogen envelope of nearly 8 M (ξ 2.5 = 7.66 × 10 −5 , ξ 1.5 = 2.38 × 10 −4 ). The H/He-layer is in the process of extended, but incomplete, 2nd dredge-up of the He-layer that is typical for AGB/sAGB stars and those at the low-mass end of the CCSN domain. The star is basically an iron-core AGB star, maybe should be called a hyper-AGB (hAGB) star.
As the envelope is not rich of metals, mass loss is expected to play only a negligible role during the star's evolution. This leaves the total mass of the star almost unchanged (see Table 1). Due to its structure it was one of the first iron-core progenitors that exploded in fully self-consistent 3D simulations by Melson et al. (2015a) with V -P and was also investigated by Radice et al. (2017) and Müller et al. (2019). The result of Melson et al. (2015a) provides the initial state for our investigation. often used also in other context, e.g., for the "α-parameter" in common-envelope (CE) studies. We also show the mass M shell contained within the corresponding radius and the post-bounce time when the outermost radius of the forward shock of our 3D models reaches the interface. E bind is the binding (i.e., internal + kinetic + potential) energy in the progenitor star outside of M map , which is the location of the final mass cut. For values see Table 3.
Moreover, we investigate a solar-metallicity M ZAMS = 9.0 M star, termed s9.0, of Sukhbold et al. (2016). Its 1.30 M iron core is surrounded by a silicon shell of 0.03 M , a carbon-oxygen layer of 0.068 M , and a helium shell of 0.169 M (see Table 1). The hydrogen envelope extends from 1.57 M up to 8.75 M (ξ 2.5 = 3.83 × 10 −5 , ξ 1.5 = 5.24 × 10 −3 ). In comparison to model z9.6, model s9.0 is just slightly less evolved on its track to 2nd dredge-up of the He core, however, the convection is better driven by the iron-peak opacity from the outside-in. The s9.0 progenitor was chosen to be representative for low-mass CCSNe by Jerkstrand et al. (2018), who studied the late-time nebular spectra of the supernova (SN), by Glas et al. (2019a) focusing on the neutrino emission during the explosion, and by Burrows et al. (2019) in a large set of 3D simulations. The three-dimensional exploding model for our investigation is provided by Melson et al. (2019) and has also been modeled with V -P . Although the progenitors considered in this study have very similar ZAMS masses, their pre-collapse core structures differ significantly (see Figure 1). We stress that the ρr 3 -profiles of the progenitors are decisive for the long-time evolution of the explosion Wongwathanarat et al. 2015). The behavior of the ρr 3 -profile yields important information on the propagation of the shock through the stellar structure, because, according to Sedov et al. (1961), positive gradients of ρr 3 cause shock deceleration, whereas negative gradients cause the opposite. Additionally, variations of the shock velocity produce crossing pressure and density gradients behind the shock front near the composition-shell interfaces (Chevalier & Klein 1978). Such conditions are essential for RT instabilities as detailed in Section 5.1, assigning them a crucial role for explaining high-velocity metal-rich ejecta and radial mixing of heavy elements (e.g., Arnett et al. 1989a;Nomoto et al. 1990;Wongwathanarat et al. 2015).
The ECSN progenitor exhibits an extremely sharp drop of the ρr 3 -profile just outside of the ONeMg-core. It is this drop in density that enables fast explosions due to an early, rapid decline of the mass accretion rate and of the corresponding ram pressure at the SN shock (Kitaura et al. 2006). Outside of the core, the ρr 3profile grows monotonically as no other composition interfaces are encountered. The z9.6 progenitor falls into the class of ECSNlike progenitors also in this respect: Similar to the electron-capture progenitor, the z9.6 model shows a monotonic growth of the ρr 3profile exterior to the CO core, where only a small step in the density profile can be seen at the He/H interface. Model s9.0 on the other hand exhibits strong variations in its ρr 3 -profile. Each interface of different composition layers is accompanied by a negative ρr 3gradient close to the interface. Of particular interest are the CO/He and He/H interfaces. These interfaces have an impact on the longtime evolution of the explosion as will be discussed in Section 5.1.
For this paper we performed spherically symmetric (1D), axisymmetric (2D) and fully three-dimensional (3D) simulations for all progenitors beyond the moment when the shock reaches the surface of the star. The different setups and approaches for these simulations will be described in the following section.

NUMERICAL METHODS AND PHYSICAL SETUP
In order to cover collapse, shock revival, and shock propagation through the envelope and circumstellar material we employ a stepwise approach. Core collapse, bounce, and post-bounce evolution until the explosion is well on its way, i.e., the neutrino-dominated phases of the first second(s) around and after core bounce, of the iron-core progenitors are simulated with the V -P code. The corresponding long-time simulations covering the shock propagation through the star, after the onset of the explosion until and beyond shock breakout, are conducted with the P -HOTB code. The explosion of the ECSN progenitor is simulated entirely with P -HOTB. 2 In the following sections we describe the numerical and physical features of the codes and the setups applied during the different stages of the evolution.

V -P code
V -P is a hydrodynamics code based on an implementation of the Piecewise Parabolic Method (PPM) of Colella & Woodward (1984), coupled with a three-flavor, energy-dependent, ray-by-ray-plus (RbR+) neutrino transport scheme that iteratively solves the neutrino energy and momentum equations with a closure determined from a tangent-ray Boltzmann solver (Rampp & Janka 2002). It employs the full set of neutrino reactions and microphysics presented in Buras et al. (2006) and the high-density equation of state (EoS) of Lattimer & Swesty (1991) with a nuclear incompressibility of K = 220 MeV. At low densities (ρ ≤ ρ HD = 10 11 g/cm 3 ) V -P uses the EoS of Janka & Müller (1995), which includes the contributions of photons, arbitrarily degenerate and arbitrarily relativistic e + /e − , and non-degenerate, non-relativistic nucleons and nuclei. The relative abundances of 23 nuclear species (including some neutron-rich nuclei) are determined by an NSE solver in regions with temperatures above T NSE . Below T NSE a flashing scheme is used to approximately treat nuclear burning (see Rampp & Janka 2002). For unshocked, collapsing stellar matter we choose T NSE = 0.5 MeV, and for neutrino-heated postshock matter we take T NSE = 0.5 MeV for the simulation of model z9.6 and T NSE = 0.34 MeV for model s9.0. 3 The simulations presented here are performed with a 1D gravitational potential including general relativistic corrections, Case A of Marek et al. (2006). The neutrino transport solver contains corrections for general relativistic redshift and time dilation effects. V -P makes use of the axis-free Yin-Yang grid (Kageyama & Sato 2004) based on the implementation of Melson (2016). 2 The reason for these different treatments is mainly technical: The V -P version used for ONeMg core collapse by Kitaura et al. (2006);Janka et al. (2008), andHüdepohl et al. (2010) has not yet been updated with the 3D developments and parallelisation optimization applied to the version used for Fe-core collapse. However, for suitable choices of parameter values, the explosion dynamics of ECSNe computed with P -HOTB is very similar to the fully self-consistent 2D V -P and C C N T-V explosion models discussed by Janka et al. (2008) and Janka et al. (2012), respectively. 3 V -P , which does not apply a nuclear network for T < T NSE . The lower value of T NSE permits us to follow the ejection of mass through neutrino heating for a longer time period in model s9.0, where the expansion velocity of the expelled matter is smaller than in z9.6, thus facilitating nucleon recombination to α particles and heavy nuclei.

P -HOTB code
P -HOTB is based on the same hydrodynamics module as V -P and uses the implementation of the Yin-Yang grid presented in Wongwathanarat et al. (2010a). It employs the EoS of Lattimer & Swesty (1991) for high densities above a threshold value ρ HD (usually 10 11 g/cm 3 ) and the "Helmholtz" EoS of Timmes & Arnett (1999) for densities below ρ HD , which takes into account arbitrarily degenerate and relativistic electrons and positrons, photons, and a set of non-degenerate, non-relativistic nuclei. The set of nuclei consists of neutrons n, protons p, 13 α-nuclei ( 4 He, 12 C, 16 O,20 Ne,24 Mg,28 Si,32 S,36 Ar,40 Ca,44 Ti,48 Cr,52 Fe,56 Ni), and an additional tracer nucleus Tr, which tracks the production of neutron-rich nuclei and replaces 56 Ni in environments with low electron fraction, Y e <0.49. These nuclear species are described as non-relativistic Boltzmann gases. The advection of the species is treated with the Consistent Multi-fluid Advection (CMA) scheme of Plewa & Müller (1999). NSE is assumed above T NSE = 9 × 10 9 K and accounted for by an NSE table including the nuclei listed above (Kifonidis 2004, private communication). Nuclear burning is considered at temperatures below T NSE with a 13-species α-network, which is consistently coupled to the hydrodynamic modeling. At the boundary between network and NSE we assume that all free neutrons and protons recombine to yield 4 He. We thus add the mass fractions of p and n onto the mass fraction of 4 He, accounting for the corresponding energy release 4 . The P -HOTB code uses a 3D gravitational potential with the general-relativistic (GR) monopole correction of Marek et al. (2006) as discussed by Arcones et al. (2007), while higher multipoles are obtained from a solution of Poisson's equation as described in Müller & Steinmetz (1995).
Different from V -P , P -HOTB uses a three-flavor grey neutrino transport scheme as presented in Scheck et al. (2006), 5 which is applicable at low and moderate optical depths. Therefore, the high-density core of the PNS, with a mass of M c = 1.1 M and densities well above those of the neutrinospheric layer, is replaced by a closed (Lagrangian) inner grid boundary at radius R ib . The excised 1.1 M PNS core is taken into account in the gravitational potential as a central point mass. The contraction of the PNS is mimicked by an inward movement of the boundary radius R ib , whose motion is followed by all grid points in the computational domain. We use the contraction of R ib as prescribed by Ertl et al. (2016). The time-dependent neutrino luminosities at R ib are imposed as boundary conditions as provided by an analytic one-zone cooling model following Ugliano et al. (2012), Sukhbold et al. (2016), Ertl et al. (2016), andErtl et al. (2020). This time-dependent treatment of the central-core region employs five parameters (p, a, R c,f , t 0 , R ib,f ; see Appendix A for definitions), which can be calibrated to yield explosions that fulfill the constraints set by observed SNe or by fully self-consistent 3D simulations of CCSNe. The reader is referred to Appendix A for a more detailed description of the parametric approach and to Table 2 for the parameters of the core model employed in our work.

Collapse and post-bounce setup in V -P
The collapse of the iron core progenitors is computed in 1D using the full set of neutrino interactions until 10 ms after bounce. Thereafter, 4 In a newer version of the code we allow only paired free neutrons and protons to recombine to 4 He, thus also satisfying charge conservation. 5 Some improvements and extensions to the neutrino transport module are provided in Appendix B. the simulations are mapped onto the three-dimensional Yin-Yang grid and random cell-to-cell density perturbations are imposed with an amplitude of 0.1%. The simulations employ a non-equidistant radial grid with initially 400 zones extending to 10 9 cm, which is refined in steps to more than 600 zones. This guarantees a resolution ∆r/r of better than 1% at the gain radius. The innermost 1.6 km are calculated in spherical symmetry to avoid time stepping constraints at the grid center. The angular resolution of the z9.6 model is 2 • . The post-bounce evolution of model s9.0 is computed with a newly implemented static mesh refinement (SMR) scheme presented in Melson et al. (2019), which increases the angular resolution to 1 • outside of the gain radius and to 0.5 • exterior to a radius of 160 km. The simulations with full neutrino transport are too expensive to continue them to late post-bounce times. At t pb 0.5 s the neutrino transport is therefore switched off and replaced by a simplified scheme for neutrino heating and cooling, which ensures an essentially seamless continuation with a minimum of transient artifacts. Details of this scheme are given in Appendix E. During this phase of simplified neutrino treatment both model z9.6 and s9.0 are simulated with uniform angular resolution of 2 • .

Collapse and post-bounce setup in P -HOTB
During the spherically symmetric simulation of the collapse up to core bounce, P -HOTB uses the parametrized deleptonization scheme described in Liebendörfer (2005). The necessary Y e (ρ)-trajectory was provided by Hüdepohl (2018, private communication) from his core-collapse simulations of the ONeMg-core progenitor e8.8 n with V -P . For the simulation of the ECSNe progenitor we take ρ HD = 10 11 g cm −3 and assume NSE in regions where the temperature exceeds T NSE = 9×10 9 K and apply the α-network for temperatures lower than this value. After bounce P -HOTB employs the grey neutrino transport scheme and modeling approach as presented in Scheck et al. (2006). Thus, the neutrino-opaque central core of the PNS is excised from the computational domain and replaced by the analytic core model of Ugliano et al. (2012). In Table 2 we list the parameter values of the PNS core model used for a set of simulations of model e8.8. We perform 1D and 2D simulations for all four sets of parameter values and choose the e8.8 10 calibration as our reference case for a 3D simulation. The 1D and 2D simulations possess a non-equidistant radial grid with 2000 zones up to a radius of R ob = 2×10 10 cm. The 3D run has only 1400 radial zones for computational efficiency. The multi-dimensional simulations are conducted with an angular resolution of 2 • , and the 3D simulation makes use of the Yin-Yang grid. We restrict ourselves to a 1D gravitational potential with GR corrections (Arcones et al. 2007) because the explosions of model e8.8 are nearly spherical.

Setup for the long-time simulations
The simulations of the long-time evolution of all models are computed with P -HOTB. For this we map the final state of a post-bounce simulation at time t map onto a new computational grid within P -HOTB, similar to the procedure described in Wongwathanarat et al. (2015). We also add the low-density extensions to the Helmholtz EoS described therein. In Table 3 we list the times of mapping and the inner and outer radii of the new computational domain, R ib and R ob , respectively. The mass contained within R ib is treated as a point mass 6 and is called M map .
The time t map is chosen such that the explosion energy has effectively converged to its asymptotic value and a neutrino-driven wind region has developed around the PNS, where the outflow is essentially spherical and reaches supersonic velocity.
All long-time simulations are computed with an angular resolution of 2 • . We use a non-equidistant (geometrically increasing) radial grid from the inner to the outer boundary. In order to guarantee sufficient resolution where needed, the radial grid is allowed to move with the ejecta starting from t pb ∼10 s. Gravity is accounted for by a 1D GR-corrected potential and nuclear reactions are still considered. When mapping the iron core models, z9.6 and e9.0, into P -HOTB, we recombine free n and p from the freeze-out of NSE into 4 He under the condition of charge conservation and account for the energy release. Moreover, we combine all neutron-rich nuclei formed in neutrino-heated ejecta into tracer (Tr) material.
When mapping from the simulations of the onset of the explosion to the follow-up simulations, the central region interior to R ib (Table 3) is removed from the computational domain. Similar to Wongwathanarat et al. (2015) we prescribe an inflow boundary condition at R ib , which corresponds to the neutrino-driven baryonic mass-loss ("neutrino-wind"; e.g Qian & Woosley 1996) generated by ongoing neutrino-energy deposition in the surface layers of the cooling PNS. In contrast to Wongwathanarat et al. (2015), we employ neutrino-wind results adopted from 1D simulations of the explosions, seamlessly connected to the fully multi-dimensional explosion simulations by choosing R ib to be in the supersonic wind region (ensuring that perturbations cannot propagate back to the inner boundary creating artifacts) and by applying the wind data at times when the outflow properties match closely between the 1D and the (angle-averaged) multi-dimensional models. This is possible because the PNSs in 1D and multi-dimensional models are extremely similar and the neutrino-emission and thus the neutrinodriven winds also have very similar properties. For the long-time run of model z9.6, we therefore employ neutrino-wind conditions of a 1D simulation of this model with the V -P code. This model treats PNS convection with an approach based on mixing-length theory and exhibits neutrino-emission properties that are hardly distinguishable from the multi-dimensional calculation (Mirizzi et al. 2016, see). The time dependence of the radial velocity v r , density ρ, and total (i.e. kinetic + internal) energy density e, which are needed for setting the boundary condition, are shown in Figure 3 (solid lines) as extracted from the 1D explosion simulation of the z9.6 model at a radius of 600 km (in the supersonic wind domain). For the 3D simulation of model e8.8 we use the neutrino-wind results from the corresponding 1D run with the same explosion energy (see dashed lines in Figure 3), whereas we do not impose a wind boundary condition in the long-time simulation of 6 We ensure that matter at radii smaller than R ib has velocities smaller than the local escape velocity and will thus eventually contribute the to final NS mass. log (e/e 0 ) Figure 3. Time-dependent behavior of neutrino-wind density ρ, radial velocity v r , and total energy density e, normalized to their initial values at t map (see Table 3), which defines the start of the long-time simulations of model e8.8 and z9.6. The data for model z9.6 are extracted from a 1D simulation of the explosion of the 9.6 M progenitor with V -P (see Mirizzi et al. 2016), evaluated at a radius of 600 km (solid lines). For model e8.8 we use data of the respective 1D simulation with P -HOTB (dashed lines). Fe + e + + ν e + γ (19%) . The above reactions provide an energy source for the surrounding plasma in the form of gamma radiation (E γ ) and kinetic energy (E kin,e + ) of the positrons that are produced in the β + decays. We include the annihilation energy of the positrons with electrons (E ann ) in E γ . The produced neutrinos escape freely. The average energy available (including the kinetic energy of the positron in the cobalt decay) per decay is E γ,Ni = 1.72 MeV, E γ,Co = 3.735 MeV (Nadyozhin 1994). A fraction of the γ's may escape depending on the (radial) optical depth τ(r) of the gas up to the stellar surface at radius R * . This optical depth is defined as where Y e the electron fraction and κ γ the optical opacity. In the practical application the integral boundary r in Equation (1) is the radial location of a considered grid cell, and we assume that the trapped fraction of the locally produced γ radiation, (1 − exp[−τ(r)]), deposits its energy locally in the same cell of the computational grid. Assuming Compton-scattering is the dominant opacity source, we adopt a constant value of κ γ = 6.0 × 10 −2 cm 2 g −1 (Swartz et al. 1995). Therefore, the energy per mass ∆E i /∆M deposited by each species i, with mass fraction X i and nuclear mass m i , into the surrounding plasma during a time step ∆t is given by where ∆X i = 1 − e −∆t/t 0, i is the change of X i during ∆t with t 0,i = t 1/2,i (ln(2)) −1 being the life-time of species i, and E kin,e + is the kinetic energy of the positron in the cobalt decay. The energy is assumed to be deposited locally, and thermodynamic quantities are self-consistently updated.

Shock propagation and explosion energetics
In the following we provide a brief overview of the most important features of the early post-bounce evolution in the 3D simulations of models e8.8, z9.6, and s9.0. The reader is referred to Melson et al. The dynamics of our simulations for the 8.8 M progenitor closely resemble these previous findings. In Figure 4 we show the evolution of the angle-averaged radius of the supernova shock and the diagnostic explosion energy of our three-dimensional simulations during the first three seconds after bounce. The angle-averaged shock radius is calculated as where dΩ= sin θdθdφ. The diagnostic explosion energy at all times is given by the integral of the total (i.e., kinetic plus internal plus gravitational) energy density in the postshock region, defined as e b = e int + e kin + e grav , over volume elements where it has a positive value The sharp drop in density outside of the ONeMg-core in model e8.8 leads to an early and strong decrease of the accretion rate and hence ram pressure at the shock. Consequently, the SN shock expands rapidly and reaches the core/envelope boundary (R He/H = 1210 km), at 0.21 s after core bounce. This is in stark contrast to the typical CCSN, where the high ram pressure stalls the shock expansion at a small radius for several 100 ms. The explosion energy starts rising steeply, fuelled by the onset of a neutrino-driven wind, as soon as the shock leaves the ONeMg-core.
The acceleration of the shock at the core/envelope boundary is followed by a drastic switch to deceleration after the shock passes the lower boundary of the H-envelope (see upper right panel of Figure 4). This is caused by a sudden change of the density gradient at the core/envelope transition. Consequently, the neutrino-heated ejecta pile up in a dense, compressed and decelerated shell behind the SN shock.
The ECSN-like structure of model z9.6 is reflected in the evolution of the SN shock front and in the growth of the diagnostic explosion energy. The shock radii remain almost perfectly spherical in both the e8.8 and z9.6 models. The acceleration of the blast wave outside of the iron core, however, ends earlier in model z9.6 than in model e8.8 because of the more gradual changes in the density profile. In both cases the explosion energies also start to rise early and saturate just after t pb ∼1 s.
In contrast to the z9.6 and e8.8 models, model s9.0 lacks the very steep density gradient at the edge of the Fe-core, as can be seen in Figure 1. The shock can expand initially up to 180 km at t pb ∼130 ms, but then it enters a phase of recession. The arrival of the Si/O interface at the shock and the decreasing mass-accretion rate within the oxygen shell eventually lead to shock expansion at ∼0.32 s after bounce. Shock expansion is aided by strong convection behind the shock front (see also Melson et al. 2019). Similar to the results presented in Glas et al. (2019a), who used the same progenitor, the model remains convection-dominated and does not exhibit any sign of the oscillatory growth of the SASI (Blondin et al. 2003;Foglizzo et al. 2007). However, although SASI does not develop in the simulation, the forward shock experiences largescale deformation with a dipole amplitude of ∼10% (compared to the angle-averaged shock radius). This deformation is driven by big plumes that form in the post-shock layer. Contrary to the ECSN-like models, strong anisotropic, non-radial mass flows persist around the PNS during several seconds after bounce. Continuous mass accretion onto the PNS through narrow funnels delays the emergence of the spherical neutrino-driven wind, which is why we needed to continue the simulations including PNS and neutrino treatment for more than three seconds after bounce. Explosion energy and shock velocity in model s9.0 remain considerably lower than in the other two progenitors (see Figure 4).

Neutrino emission properties
In the following we present the neutrino emission properties of our two 3D simulations performed with the V -P code. Figure 5 displays the neutrino luminosities (defined as 4πintegrated energy fluxes) and mean neutrino energies (defined as ratio of angle-averaged energy density to number density) for ν e ,ν e and heavy-lepton neutrinos ν x , as well as the three lowest-order multipoles (monopole, dipole, and quadrupole) of the electron-neutrino lepton-number flux for models z9.6 and s9.0 as functions of postbounce time (z9.6 left column, s9.0 right column). All quantities are evaluated at 400 km (transformed to an observer frame at rest at infinity). The formulas for the spherical harmonics decomposition are provided in Appendix C.
In the case of z9.6 the luminosities of all three species become very similar after only ∼180 ms, signalling the end of PNS accretion caused by the quick onset of the explosion and the rapid shock expansion. In contrast, PNS accretion continues at a significant rate until roughly 350 ms in model s9.0. Only afterwards the luminosities in this model converge to nearly the same level, mirroring the characteristic trend when the cooling emission of the newly formed NS begins. Overall, the time evolution and the values of the neutrino luminosities and mean energies of both models are very similar to each other, consistent with the nearly equal masses of the PNSs in both cases (see Table 4). Model s9.0 exhibits additional accretion emission between ∼150 ms and ∼350 ms, which enhances the ν e andν e luminosities slightly and drives a continuous increase of the mean energies of ν e andν e . However, the trend does not persist for long enough (and the PNS does not gain enough mass) to reach a crossing of the electron antineutrino energy with the heavy-lepton neutrino energy as reported by Marek et al. (2009) for more massive progenitors with more massive PNSs.
The bottom panels of Figure 5 demonstrate that both models develop a lepton-number emission dipole that dominates the higherorder multipoles and can reach and even exceed the monopole for longer periods of time (i.e., hundreds of milliseconds, see also Tamborra et al. 2014b). Note that the dipole amplitude displayed in Figure 5 is, by normalization, one third of the dipole amplitude considered by Tamborra et al. (2014b) (see also Appendix C). A long-lasting lepton-number emission dipole, whose direction is nearly stable or migrates only very slowly, is a characteristic feature of the LESA (Lepton-Emission Self-sustained Asymmetry) phenomenon that was first witnessed in 3D V -P simulations with neutrino transport by Tamborra et al. (2014b) and Janka et al. (2016). This striking phenomenon has meanwhile been confirmed with fully multi-dimensional instead of ray-by-ray neutrino transport by O'Connor & Couch (2018b), Glas et al. (2019b) and Vartanyan et al. (2019).
Neither in model z9.6 nor in s9.0 does SASI play a role and, in addition, model z9.6 explodes after a very brief period of postbounce accretion. In both models the PNS emission is therefore not masked by asymmetric accretion due to SASI shock sloshing or spiral motions. Because accretion in particular in model z9.6 does not contribute to the neutrino emission at any significant level after the onset of the explosion, its lepton-number flux asymmetry is merely determined by asymmetric convection inside of the PNS, where outward transport of lepton number is strongly suppressed in the anti-LESA direction. Therefore the lepton-number flux can even be negative in the anti-LESA direction (see Fig. 6, two upper left panels). This means that one hemisphere of the PNS radiates a greater number ofν e than ν e , whereas there is the usual excess of ν e number loss from the other hemisphere. The models, in particular z9.6 with no long-lasting accretion, therefore confirm that LESA is a phenomenon primarily generated by hemispherically asymmetric convection in the interior of the PNS, well below the neutrinosphere; the reader is referred to the more detailed discussion by Glas et al. (2019b), where also a 3D simulation with successful explosion of  We also show the total PNS kick velocities (thick lines) and the hydrodynamically induced PNS kick velocities (thin lines) in the lower right panel. For better visibility, the inset of this panel displays the PNS kick velocities caused by asymmetric neutrino emission due to the LESA phenomenon (see text for details). The total velocities are the vector sum of the hydrodynamic PNS kick and the neutrino-induced kick. Despite nearly equal neutrino-induced kicks in z9.6 and s9.0 and lower hydrodynamic kick in z9.6, this latter model has a higher total kick velocity for some time, because its hydrodynamic and neutrino-induced kick directions are essentially parallel, in contrast to the situation in model s9.0. For the iron-core progenitors we can track only hydrodynamic contributions to the kick after the transport module of V -P is switched off at t neut = 0.45 s and t neut = 0.49 s for models z9.6 and s9.0, respectively. For the ONeMg case our simplified neutrino treatment with the excised core of the PNS does not allow us to monitor the LESA induced kick. The kick of model e8.8 is scaled by a factor of 5 for better visibility. the 9.0 M progenitor is evaluated. 7 It is also interesting to note that the lepton-emission multipole ( 0) that grows fastest and thus develops high amplitudes of asymmetry first is the one of order = 4, which is also in line with the results reported by Glas et al. (2019b).
Interestingly, there are phases in both of our models when the dipole and quadrupole amplitudes become similar ( Figure 5). This is also suggested by the Aitoff projections for time t pb =0.45 s near the end of our simulation of z9.6 ( Fig. 6). At this late time a quadrupole pattern is superimposed on a hemispheric dipole asymmetry, whereas at t pb = 0.30 s a much cleaner dipole is present. This reflects the evolution of the multipole amplitudes visible for z9.6 in the bottom left panel of Figure 5. The late-time rapid growth of the dipole in the lepton-number emission of model s9.0 at t pb > 0.43 s (bottom right panel of Figure 5) is an apparent feature caused by a transient phase of a reduced lepton-emission dipole between ∼0.35 s and ∼0.45 s. This reduction is a consequence of an accretion asymmetry that channels matter towards the PNS predominantly in one hemisphere, whereas a mighty outflow in the form of a huge, rising bubble develops in the opposite hemisphere (see the bottom panels of Figures 7 and 8). The main downflow direction is misaligned with the LESA dipole direction by roughly 90 • and fuels enhanced lepton-number emission by the one-sided accretion. This enhanced lepton-number emission combined with the displaced LESA dipole strengthens the quadrupole component of the electron lepton-number emission and at the same time weakens its dipole. The remaining dipole level in the interval of [0.35 s, 0.45 s] signifies that the LESA dipole is stronger than the emission of lepton number by accretion. The dipole amplitude recovers to the full LESA Table 4. Overview of PNS properties in our multi-dimensional models at t map . Notes: All values are given at the end of our explosion simulations with neutrino treatment (t map ). E exp is the diagnostic explosion energy, which is essentially identical to the final explosion energy because of the small envelope binding energy (see Table 5). v tot NS is the total NS kick velocity resulting from the hydrodynamic (v hyd NS ) plus the neutrino-induced (v ν NS ) contributions. We measure the contribution of neutrinos to the total kick until the neutrino transport is switched of at t neut = 0.45 s and t neut = 0.49 s after core bounce for models z9.6 and s9.0, respectively. Hydrodynamic contributions and total kick velocities are monitored until t map . J NS is the total angular momentum transported to the PNS through a radius of 100 km until t map . θ v J is the angle between the direction of the total kick velocity v v v and the direction of J J J. α ej and α ν are the final hydrodynamic and neutrino anisotropy-parameters, respectively (see Appendix D). For the iron-core progenitors we average α ν over the time when the emission dipole is largest until the end of the simulation. M b is the baryonic PNS mass, which is defined as the enclosed mass within the radius R NS , at which the density drops below 10 11 g cm −3 . Note that for the simulations of the e8.8 progenitor, R NS is determined by the chosen parameters of the inner grid boundary (see Table 2). M map (see also Table 3) is the central mass contained within the excised region, from which we calculate the gravitational mass M g for a PNS radius of 12 km (see Appendix D). P NS is the PNS spin period at the end of our post-bounce simulations, assuming a final NS radius R NS of 12 km, angular momentum conservation, and a gravitational mass of M g . strength after t pb ∼0.45 s because the mass-accretion rate onto the PNS declines continuously.
The Aitoff projections show that not only the electronneutrino lepton-number flux exhibits large-scale (low-order multipolar) asymmetries, but also the ν e plusν e energy flux as well as the total neutrino energy flux (the summed energy fluxes of ν e plus ν e plus four times that of ν x ). However, while the directional variations of the lepton-number flux can be several times bigger than the angular average of this quantity, the ν e plusν e energy flux varies only within roughly 6-8% and the total neutrino energy flux even less within about 4-6% ( Figure 6). These directional variations, in particular the hemispheric asymmetries, have consequences for NS kicks and the neutron-to-proton ratio in the ejecta. This will be discussed in the following Sections 4.3 and 4.4.

Neutron star kicks by asymmetric mass ejection and neutrino emission
In order to compare the final state of the post-bounce simulations of our model set, we present in Figure 7 planar slices showing the entropy at the times t map when we start our long-time simulations. The left panels are chosen such that they align with the plane of the largest shock deformation. The right panels align with the plane of the smallest shock deformation. Note the basically spherical shape of the shock and the mild post-shock asymmetries in both the e8.8 and z9.6 models. In strong contrast to the ECSN-like models, model s9.0 shows a clear dipolar shock deformation and ejecta morphology.
The asymmetries that develop during the explosion also affect the kick of the PNS. In Table 4 we provide properties of the PNS resulting from our post-bounce simulations. 8 We list the PNS kick velocity (v NS ), ejecta anisotropy parameter (α ej ), PNS angular momentum (J NS ), angle between total PNS kick and angular momentum vector θ vJ and PNS spin period (P NS ) along with the PNS radius (R NS ) and baryonic PNS mass (M b ) at the moment in time (t map ) when we terminate the post-bounce simulations that include neutrino transport or neutrino heating (see also Table 3). The PNS radius R NS is defined as the radius where the angle-averaged density drops below 10 11 g cm −3 . The PNS mass M b is the mass contained interior to R NS . M map is the mass contained within our inner grid boundary, and thus removed from the hydrodynamic grid at the start of the long-time simulations at t map (see Table 3). M g the corresponding gravitational mass.
The acceleration of the PNS is caused by two different mechanisms. Firstly, aspherical ejection of matter leads to a gravitational tug, which accelerates the PNS into the hemisphere opposite to the maximum shock expansion and fastest ejecta Wongwathanarat et al. 2010bWongwathanarat et al. , 2013. Secondly, anisotropic neutrino emission, due to the LESA phenomenon (Tamborra et al. 2014b), can accelerate the PNS opposite to the direction of the largest total neutrino-energy flux. LESA manifests itself in a dominant and stable = 1 spherical harmonics mode of the lepton-number emission and a corresponding energy-emission dipole amplitude of several percent compared to the monopole (see Tamborra et al. 2014a;Tamborra et al. 2014b, and Section 4.2). LESA is observed in both simulations conducted with V -P . The almost spherical explosions of the ECSN-like progenitor yield very low hydrodynamic kick velocities by the "gravitational tug-boat effect" (Gessner & Janka 2018). Anisotropic neutrino emission cannot be evaluated in our simulation of model e8.8, because of the spherical treatment of the central PNS region. The kick contribution by anisotropic neutrino emission in model s9.0 is of a magnitude comparable to the contribution associated with the aspherical ejection of matter and even exceeds the hydrodynamic kick contribution in model z9.6.
We show in Figure 8 planar slices of the electron fraction (left panels) and entropy per nucleon (right panels) of the iron-core progenitors. The black, blue, and green arrows indicate the total, hydrodynamic, and neutrino-induced directions of the PNS kick.
In the case of model z9.6, the hydrodynamic and neutrino-induced kicks are nearly aligned. This results from the anti-correlation of the LESA lepton-number emission dipole and the direction of the maximum ν e +ν e as well as total neutrino luminosity (see Tamborra et . Neutrino luminosities (top), mean energies (middle), all averaged over all viewing directions, and lowest-order multipole moments (A 0 , A 1 , A 2 for monopole, dipole, and quadrupole) of the electron-neutrino lepton-number flux (bottom) as functions of time for models z9.6 (left) and s9.0 (right). All quantities are evaluated at 400 km, transformed to an observer frame at rest at infinity. With ν x we denote one species of the heavy-lepton neutrinos. Note that the dipole moment, A 1 , plotted here is one third of the dipole amplitude of the lepton-number flux defined by Tamborra et al. (2014a) (see Appendix C for details).
2014b and Section 4.2). Increased heating by ν e +ν e absorption in the postshock region in the hemisphere opposite to the LESA dipole pushes the supernova shock to larger radii. The induced asymmetry of the post-shock ejecta leads to the hydrodynamic acceleration of the PNS in the opposite direction due to the gravitational pull of the slower ejecta. The neutrino-induced NS kick acts in the same direction, i.e., nearly aligned with the LESA dipole, because of the maximum total neutrino luminosity being in the anti-LESA direction.
Due to the weak postshock convection and only small-scale ejecta asymmetries in model z9.6, the anisotropic neutrino emission produces the dominant contribution to the PNS kick. While the gravitational tug of the ejecta can only account for v hyd NS ∼ 10 km s −1 , the asymmetric emission of neutrinos alone would cause a kick velocity of v ν NS ∼ 25 km s −1 . The total kick sums up to v tot NS ≈ 35 km s −1 . A correlation or even alignment of the neutrino-induced and hydrodynamic kicks is less clear in model s9.0. The asymmetric neutrino energy emission is also responsible for a sizeable contribution to the PNS kick in this case, and even for the dominant contribution during the first 0.5 s after bounce. However, the neutrino-induced kick in s9.0, despite being in the same hemisphere as the LESA dipole vector, is not closely aligned with the LESA dipole direction as in model z9.6. The hydrodynamic PNS kick, however, is in the opposite hemisphere, different from the case of model z9.6. Both, this different orientation of the hydrodynamic kick relative to the LESA dipole, and the misalignment of LESA direction and neutrino-induced kick in model s9.0, are a consequence of the fact that hydrodynamic instabilities in the postshock layer are much Figure 6. Full-sphere Aitoff projections of various LESA related quantities for model z9.6 at 0.30 s and 0.45 s after bounce. The two upper left panels display the relative variation of the electron-neutrino lepton-number flux, (F n νe − F n νe )/ F n νe − F n νe (normalized by the angular average) at a fixed radius of 400 km as function of polar angles θ and φ. The two upper right panels show the variation of the electron fraction Y e in the neutrino-heated ejecta at 250 km. The two lower left panels visualize the relative variation of the energy flux of ν e plusν e , and the two lower right panels the corresponding relative variation of the total energy flux as sum of all contributions of ν e ,ν e and ν x , again evaluated for lab-frame quantities at 400 km. Note that the amplitude of the energy-flux variation is considerably lower than that of the lepton-number flux, and the dipole directions of both fluxes possess opposite orientations. stronger and longer-lasting. Therefore the corresponding asymmetries of the mass distribution around the PNS are more extreme in model s9.0 than in z9.6. The direction of the neutrino-induced kick is thus affected by the neutrino-emission dipole that is associated with one-sided PNS accretion. This asymmetric neutrino-emission due to accretion is transiently superimposed on the LESA dipole during the time interval of [0.35 s,0.45 s] after bounce. It is displaced by roughly 90 • from the LESA direction (see Section 4.2) and thus shifts the neutrino-induced kick away from the LESA direction ( Figure 8, lower panels).
The neutrino-heating asymmetry associated with LESA, causing more heating in the anti-LESA direction (where the ν e plusν e luminosity is higher and also the number flux of theν e with their harder spectra is higher), has a mild influence only in model z9.6. Therefore the explosion in this model is slightly stronger in the anti-LESA direction, and the hydrodynamic PNS kick is nearly aligned with the LESA vector and with the direction of the neutrino-induced PNS kick. In contrast, in model s9.0 the neutrino-heating asymmetry associated with LESA is not powerful enough to determine the deformation of the explosion. Convective mass motions, the corresponding accretion asymmetries, and the associated anisotropic neutrino emission and absorption between PNS and SN shock play a more important role and have a dominant influence on the morphology of the postshock flow. In model s9.0 the shock expansion and explosion are weaker in the anti-LESA direction, where convective downdrafts towards the PNS are numerous and massive and thus direct the hydrodynamic PNS kick towards this hemisphere (blue arrows in the southern hemisphere of the lower panels of Figure 8). As explained above, accretion-associated asymmetries of the total neutrino luminosity are also the reason why the neutrinoinduced kick is not well aligned with the LESA lepton-emission dipole. The combined hydrodynamic and neutrino-induced kicks in model s9.0 lead to a PNS velocity of ∼41 km s −1 at mapping time t map . It must be emphasized that in both models, z9.6 and s9.0, the V -P neutrino transport was switched off after about 0.5 s of post-bounce evolution, namely at ∼0.45 s and ∼0.49 s after bounce, respectively. At this time the neutrino-induced kicks have not reached their final values, which might be higher than the numbers presented in Table 4.

Electron fraction of neutrino-heated ejecta
Since the explosions of models z9.6 and s9.0 are computed with detailed neutrino transport, we describe the evolution and distribution of their electron fraction, Y e , in the following. Corresponding 2D results obtained for ECSN models in simulations with the V -P and C C N T-V codes were presented by Wanajo et al. (2011) and Wanajo et al. (2018), respectively. The electron fraction in the neutrino-heated ejecta of our 3D model e8.8 is only approximate because of the simplified treatment of the neutrino physics and the replacement of the central 1.1 M PNS core by an inner grid boundary in the simulations with the P -HOTB code (see Section 3.2). Therefore we will not further discuss the neutron-to-proton ratio in the ejecta of this model but refer the reader to Wanajo et al. (2011) and Wanajo et al. (2018).
The Y e distributions of models z9.6 and s9.0 in the left panels of Figure 8 exhibit a clear asymmetry: proton-rich neutrino-driven outflow (red and orange) in the hemisphere of the LESA dipole vector, and neutron-rich neutrino-wind ejecta (blue hues) in the opposite hemisphere. The proton excess is a consequence of the fact that the PNS develops an enhanced emission of ν e relative toν e in the hemisphere of the LESA direction, which leads to a predominant production of protons by ν e absorption on neutrons (ν e + n → e − + p) in the ejecta of the same hemisphere. In the anti-LESA direction the emission ofν e by the PNS is relatively higher, allowing for more efficient creation of neutrons through the reactionν e + p → e + + n in the neutrino-heated outflow.
In addition to this spatial asymmetry in the neutron-to-proton ratio, which is present in both explosion models computed with the V -P code, model z9.6 exhibits neutron-richness in the mushroom-shaped heads of Rayleigh-Taylor plumes in all directions, very similar to previous findings in 2D explosion simulations of the low-mass progenitors of models e8.8 and z9.6 with self-consistent neutrino transport (see Wanajo et al. 2011;Wanajo et al. 2018). The plumes contain the earliest neutrino-heated ejecta. Their neutron excess is explained by the extremely fast outward propagation of the SN shock in these two models, which is enabled when the mass-accretion rate of the shock plummets because of the steep density decline at the edge of the degenerate core. The fast shock acceleration allows the neutrino-heated, buoyant gas to expand very quickly away from the gain radius. The gas rises so rapidly that neutrino absorption is unable to increase Y e from its low values (∼0.25-0.35) near the gain radius to 0.5 or higher. Instead, the mushroom plumes stay considerably neutron rich in model z9.6, in contrast to model s9.0, where the slower shock acceleration and the correspondingly slower expansion of the first bubbles of neutrino-processed matter lead to Y e around 0.5 due to the dominant absorption of ν e (Figure 8).  Figure 9. Left panels: Distribution of Y e in cross-sectional slices of 2D (left) and 3D (right) simulations of model z9.6 at 0.13 s (top) and 0.4 s (bottom) after core bounce. The ejection of neutron-rich matter (with 0.39 Y e 0.50) in fast-rising, buoyant plumes is visible. Note that the chosen color range (lower color bar) does not allow to display the Y e variation at radii below ∼100 km. The colored lines visualize the infall and escape trajectories (projected onto the x-y-plane of the cross-sectional slice, the flow direction is marked by arrows) for three selected mass elements (trajectories A,B,C), which ultimately get expelled in the neutron-rich mushroom heads. The local Y e value along the trajectories is color coded according to the upper color bar. The Y e evolution of these fluid elements is more quantitatively represented in the right panels. The locations of these elements at the time of the plot (0.13 s) are indicated by small black circles. The white dashed circle marks the mean neutrinosphere, defined by a spectrally averaged optical depth of 2/3. The lower left panels, besides showing the farther expanded, neutron-rich plumes, display neutrino-driven wind ejecta, which develop neutron-rich conditions in one hemisphere (bluish colors) and proton-richness in the opposite hemisphere (reddish colors) because of the LESA dipole asymmetry of the ν e andν e emission from the nascent neutron star. We stress that the increase of Y e in the expanding mushroom heads between 0.13 s and 0.40 s is not a consequence of neutrino reactions (which basically cease at r 250 km). Instead, a gradual, slow Y e increase is caused by numerical diffusion connected with the ejecta mass flowing over the Eulerian grid. Upper right panel: Evolution of the electron fraction for the three selected mass elements that are ejected in the neutron-rich mushroom heads. The Y e evolution during the infall is shown by thin solid lines, during the re-ejection by thick solid lines. Note that during the infall all three lines, A, B, and C, lie on top of each other, starting at Y e = 0.5. The dashed lines represent the local electron fraction, Y equil e , for reactive equilibrium. Bottom right panel: Expansion timescale, τ exp (solid lines), and timescale of Y e changes, τ dYe /dt (dashed lines), measured along the outflow parts of the solid curves in the upper right panel. plane) of three infalling, then neutrino-heated, and finally outgoing mass elements. The color-coding of these trajectories as well as the upper right panel show that the infalling matter starts with a value of Y e = 0.5, gets neutronized by electron captures during infall, reverses direction at a minimum radius around 100 km close to the gain radius, and experiences an increase of Y e again as it speeds away from the gain radius. The rise of Y e , however, flattens when the expansion time scale, τ exp = r |v r | −1 (with v r being the local radial component of the velocity), becomes shorter than the timescale of Y e changes, τ dY e /dt = Y e Q Y e −1 (with Q Y e being the local net rate of ν e andν e emission and absorption reactions as computed during the hydrodynamic simulation). This can be seen in the upper and lower right panels. The electron fraction in the expanding ejecta stays well below the value for local kinetic equilibrium, Y equil e , which corresponds to the condition when ν e andν e absorption and emission reactions are balanced for the density and temperature of the tracked mass elements at each radius r. Close to the point of return the radial velocities become small and τ exp increases steeply, whereas farther out |v r | ∝ r roughly holds and τ exp becomes nearly constant. The values of τ dY e /dt for trajectories A and C increase at small radii, because these mass elements approach the NS so closely that their Y e at the return point drops to values near Y equil e , for which reason Q Y e , the net rate of Y e changes, becomes low. In contrast, mass element B returns outward at a considerably larger radius and its Y e value never gets close to the local Y , the corresponding panels above and to the right show the marginal distributions. Blue and red bars indicate Y e < 0.5 and Y e > 0.5, respectively. In order to minimize effects of numerical diffusion, the mass distributions for Y e and s were measured for outflowing material (v r > 0) at a radius of r = 250 km.
expansion timescale, because all neutrino rates become slow and ultimately cease when the mass elements gain distance from the neutrino source and move to low densities and temperatures. This explains why the asymptotic values of Y e along the trajectories stay well below the local equilibrium values at large radii.
In order to unravel similarities and differences between our current 3D models and the 2D explosion simulations of low-mass progenitors considered previously by Wanajo et al. (2018), we compare the ejecta properties of our 3D models to corresponding 2D results in the following three paragraphs. Neutron-rich, buoyant, fast plumes are not only present in the 3D case but also in the corresponding 2D models (left panels in Figure 9). However, in the 2D case extended regions around the PNS with proton excess or neutron excess exist in all directions in the northern as well as southern hemisphere. If this finding is connected to a 2D phenomenon that corresponds to LESA in 3D, the LESA direction is less stable in 2D than in 3D. Indeed, an inspection of the dipole of the electronneutrino lepton-number flux shows that the dipole direction in 2D flips from one hemisphere to the other with a full cycle period of ∼0.1 s. Moreover, the grid axis seems to have a disturbing influence that leads to artificial effects, because both in the northern and southern directions proton-rich, collimated outflows appear very close to the axis, where they are much stronger than at large angles away from the axis. Because of the influence of the artificial Table 5. Total masses of chemical elements M i of the pre-collapse progenitor models as shown in Figure 2 (top), behind the shock front in our 3D SN simulations at the time of mapping to the long-time runs, t map , and at the time when the entire shock has broken out from the surface of the exploding star, t sbo . Also given are the binding energies still ahead of the shock at time t map .
Notes: n and p label free neutrons and protons, respectively, left unbound after the freeze-out from NSE, whereas H labels hydrogen from the envelope. Since the onset of the explosion of model e8.8 was simulated with P -HOTB, no free protons remain in the ejecta. The high mass of 4 He compared to 56 Ni and Tr in model z9.6 at t map is a consequence of the high value of the temperature when NSE was switched off in this simulation, T NSE ∼ 5.8×10 9 K, instead of T NSE ∼ 4.0×10 9 K in model s9.0. E bind (r > R sh ) is the total binding (i.e. internal + kinetic + gravitational) energy ahead of the shock at the mapping time t map . Free neutrons at t map are assumed to decay and are added to free protons at t sbo . Since our network in the long-time runs contains only 56 Ni, 56 Co, and 56 Fe, all other iron-group species (from the explosion models or the progenitor) are included in the mass of these nuclei listed for t sbo . symmetry axis, it is not easy to diagnose whether the flipping 2D dipole has any physical relation to the stable or slowly migrating lepton-emission dipole of LESA in 3D.
In Figure 10, we present the ejecta mass distributions versus entropy per nucleon, s, and electron fraction, Y e , for 2D and 3D simulations of z9.6 (top) and s9.0 (bottom), including the fast convective plumes as well as the neutrino-driven wind material. The mass distributions are constructed by integrating all matter that flows through a sphere of 250 km radius with positive radial velocities. This choice of radius ensures that neutrino interactions have essentially ceased at this location, but effects due to numerical diffusion and mixing in the outflowing material, when it moves across the Eulerian (i.e., spatially fixed) computational grid, are minimized in the mass-versus-Y e distributions extracted from the simulations for nucleosynthesis discussions. A comparison of the upper and lower left panels in Figure 9 demonstrates the consequences of this diffusion. The mushroom heads are very neutron rich at 200-300 km and t pb = 0.13 s. However, although the neutrino reactions basically cease at r 250 km, the electron fraction in the nearly self-similarly expanding mushroom heads still continues to increase slowly to values closer to 0.5 (see lower left panels for t pb = 0.4 s. This is partly an unphysical, but unavoidable, consequence of the fact that the mushroom heads experience gradual numerical mixing (besides some unclear degree of physical mixing) with the surrounding postshock matter, which possesses Y e = 0.5.
The 2D and 3D mass distributions exhibit close similarity in their shapes as well as widths, with differences only in details of their substructure. These differences originate mainly from the fact that in 2D there are only a few plumes, which are axially symmetric, massive objects, whereas in 3D there is a greater number of such bubbles, of which each one contains less mass than the toroidal 2D objects. The 3D distributions are therefore smoother and possess less fine structure. The total masses of neutrino-processed ejecta until about 0.5 s after bounce are 10.34×10 −3 M in model z9.6 and 2.64×10 −3 M in model s9.0. In both cases about 60% of these masses are in the fastest mushroom-shaped plumes, the remaining 40% are carried by the neutrino-heated subsequent outflows. Also in both cases roughly 50% of this latter component are neutronrich (1.90×10 −3 M in z9.6 and 0.55×10 −3 M in s9.0), and about 50% are proton-rich (2.31×10 −3 M in z9.6 and 0.51×10 −3 M in s9.0).
The mass-versus-Y e distributions of model z9.6 in 2D and 3D resemble closely the corresponding distributions obtained in 2D explosion simulations for the ONeMg-core progenitor e8.8 n (mentioned in Section 2.1) and for our 9.6 M Fe-core progenitor and another ultra-metal-poor Fe-core progenitor of 8.1 M , whose results were presented in Wanajo et al. (2011) and Wanajo et al. (2018) (see figure 2 in both papers). In all cases the distributions are very wide, stretching from Y e ∼ 0.55-0.6 on the proton-rich side to Y e ∼ 0.40 (or even a bit lower in the 2D models) on the neutron-rich side. The corresponding nucleosynthetic abundance patterns are extremely similar and characteristic of this class of "ECSN-like" explosions, with high production factors for light trans-Fe elements from Zn to Zr, resulting from the appreciable ejection of neutronrich matter in these models (Wanajo et al. 2018). In contrast, model s9.0 displays a mass-vs-Y e distribution that resembles the result for the 11.2 M model s11 in Wanajo et al. (2018). It is more strongly peaked around 0.5 and has steeper and less extended left and right wings. In particular, considerably less matter with neutron excess is ejected, and the minimal Y e is around 0.46 instead of 0.40. The reason for this difference are the lower velocities with which the neutrino-heated ejecta expand outward from the gain radius. This slower expansion allows ν e -absorption to lift Y e from its low initial values at the gain radius to values closer to 0.5. Such conditions are less favorable for the production of neutron-rich trans-Fe species. The abundance pattern in the neutrino-processed ejecta of model s9.0 must be expected to show similarities with model s11 of Wanajo et al. (2018). A detailed investigation of the formation of chemical elements in s9.0 in comparison to z9.6 is deferred to future work.

Elemental distribution shortly after shock revival
In order to assess the strength of physical mixing due to multidimensional hydrodynamical flows during the first seconds in our 3D models, we present in Figure 11 normalized mass distributions for selected chemical elements behind the shock radius as functions of velocity (top panels) and mass coordinate (lower panels). To sample the velocity space, we choose 50 bins between the maximal and minimal velocities within the considered region. We use 30 bins to sample the distribution in the mass coordinate. The total masses in the postshock volume at the time t map when we map our models onto the new computational grid for the long-time simulations are listed in Table 5.
The quasi-spherical explosion of model e8.8 manifests itself also in the distribution of the chemical elements over mass and velocity coordinates, since the initial shell structure of the progenitor is essentially preserved. Because of the extreme shock acceleration in the steep density gradient at the edge of the degenerate core, most material of the thin carbon shell of the progenitor travels with 30, 000−70, 000 km s −1 ahead of the outer mass shells of the former ONeMg-core, which propagate with up to ∼60, 000 km s −1 . Most of the newly synthesized iron-group and α-nuclei expand with velocities below ∼30, 000 km s −1 . Inefficient mixing also confines the neutrino-heated ejecta within M(r) ≤ 1.338 M .
Slightly more powerful convective overturn in model z9.6 leads to more efficient mixing in mass and velocity space (see middle panels of Figure 11). Therefore some neutrino-heated material gets mixed into the carbon and oxygen shells and travels at more than 20, 000 km s −1 . The bulk of the metal-rich ejecta still expands with slower velocities, however. Similar to model e8.8, the deceleration of the SN shock at the CO/He interface in model z9.6 also compresses parts of the ejecta into a dense shell. Model s9.0, again, differs distinctively from the ECSN-like models. Strong and longlasting convective overturn in the postshock layer completely erases the initial onion-shell structure of the progenitor. The chemical elements become nearly homogenously mixed over mass and velocity coordinates (see right panels of Figure 11). The fastest neutrinoheated ejecta are associated with a big high-entropy plume inducing a dipolar deformation of the shock wave, which can be seen in Figure 7.

Linear stability analysis
As mentioned already, the velocity of the forward shock depends on the progenitor structure, in particular the ρr 3 -profile. According to Sedov et al. (1961) one expects an increase/decrease in the shock velocity when the gradient of the ρr 3 -profile is negative/positive. This acceleration and deceleration lead to density and pressure gradients in the post-shock region with opposite signs. Thus, perturbations in the matter distribution become unstable to the Rayleigh-Taylor (RT) instability (Rayleigh 1882;Chevalier & Klein 1978). As the ρr 3profiles, and hence the shock speed, vary significantly between the progenitors, the growth of RT instabilities will affect the long-time evolution of our models in different ways.
In order to aid us with the interpretation of our threedimensional simulations, we appeal to 1D long-time simulations which are started from angle-averaged states of the 3D explosion models at the same times as given in Table 3. The numerical and physical setup for the spherically symmetric simulations remains unchanged compared to the 3D runs. The 1D models can teach us about the behavior of the shock, while it propagates through the envelope. Additionally, we can compute the linear RT growth rates σ RT of small initial perturbations by tracking Lagrangian mass coordinates in our 1D simulations, following Müller et al. (1991a). In the incompressible case the growth rate is given by where p and ρ are the pressure and density. In the compressible case the growth rate becomes where c s is the speed of sound and γ the adiabatic index. Equation (6) is less restrictive than its incompressible counterpart 9 . For the time-dependent amplification factor of an initial perturbation ξ 0 we integrate Equation (6) according to This analysis enables us to estimate the time and locations in mass coordinate where the fluid becomes unstable to the RT instability, helping us to understand the origin of outward/inward mixing of different layers of chemical elements. In Figure 12 we show the amplification factors in the incompressible (colored lines) and compressible (gray line) cases for our spherically symmetric simulations. Significant differences between the models are evident.
The left panel of Figure 12 shows the integrated growth rate of the 1D version of model e8.8 10 . The extreme deceleration of the forward shock when moving into the hydrogen envelope induces very high growth factors at M(r) ∼ 1.34 M . The ECSN-like structure of model z9.6 is reflected in its amplification factors. Strong deceleration of the forward shock at the CO/He interface creates the necessary condition for the RT instability to grow there. As can be seen in Figure 12, the amplification factors at ∼1.37 M reach extremely high values as observed in the ECSN case, too. Interestingly, no growth is expected at the He/H interface, which is a consequence of the tiny step in the ρr 3 profile at this interface (see Figure 1, bottom left panel, blue curve, at r ≈ 10 12 cm).
Model s9.0 shows striking differences to the models described 9 Note that we are interested in the locations of maximal growth but not in the actual values. Calculating the growth factors from multidimensional models using angle averaged-values gives similar overall structure but lower amplitudes  s9.0 3.14 s Figure 11. Normalized binned distributions of chemical elements as functions of radial velocity (top panels) and enclosed mass (bottom panels) for all of the unbound postshock material in our 3D models at t map (see Table 3). We use 50 bins in velocity space and 30 bins in the mass coordinate. The mass coordinate in the 3D models is defined by the mass enclosed by given radius and starts at the mass value contained by the PNS.
above. Due to its more shallow ρr 3 -profile in the core region and the consequently weaker episodes of shock deceleration and acceleration, the peak amplitudes of the growth factors are smaller. They are, however, of similar magnitude as in the more massive models investigated by Wongwathanarat et al. (2015). Two distinct regions of instability can be discerned around the CO/He and He/H interfaces. Similar to the RSGs presented in Wongwathanarat et al. (2015), the strongest contribution arises from the He/H interface followed by the CO/He interface, where the amplification factors are about 10 orders of magnitude lower. While the compressible and incompressible analyses give similar results for models e8.8 and z9.6, the compressible evaluation of the growth factors in model s9.0 predicts additional growth within the He-shell of the star. This is caused by the passage of the reverse shock at later times (t pb > 3×10 3 s). Note, however, that at this advanced stage of evolution the instability is already in a strongly non-linear regime, where Equation (6) loses its validity.

Propagation of the supernova shock
In Figure 13 we show the angle-averaged shock radii R sh , reverseshock radii R rsh and shock velocities v sh for our 3D long-time simulations.

3D model e8.8
In model e8.8 the forward shock has passed the He/H interface and is traveling through the H-envelope. It left the core with about 100, 000 km s −1 and is progressively decelerated due to the monotonically increasing ρr 3 within the H-envelope (see also Figure 4). Since the innermost neutrino-heated ejecta are not in sonic contact with the SN shock, they do not feel the deceleration and thus travel with constant velocity of ≈35, 000 km s −1 until they catch up with the immediate post-shock matter at t pb ≈ 130 s. The interaction of the fast neutrino-heated ejecta with the post-shock material keeps the shock at slightly higher velocities (see small kink at r ∼ 10 7 km, lower left panel in Figure 13). Nevertheless, the forward shock decelerates continuously (without relevant phases of acceleration) until it reaches the surface of the star. Note the dramatic deceleration of the shock from its initial velocity of 75, 000 km s −1 at t pb = 2.55 s down to <1000 km s −1 when it leaves the star. A large fraction (∼75%) of the kinetic energy of the shock wave is used up for heating the hydrogen in the outer layers of the envelope, leaving only around 2.5×10 49 erg of kinetic energy in the ejecta.
The early strong deceleration of the forward shock at the bottom of the H-envelope at t pb 0.25 s (see also Figure 4) leads to the compression of the CO and He layers into a high density shell. Less than 1 s later, a strong reverse shock forms at the base of this shell as the forward shock is progressively slowed down. The dense shell separates the postshock material from the neutrino-heated wind ejecta, which can be seen in the density profiles shown in Figure 14.
It is remarkable that in the ECSN the neutrino-heated ejecta are denser than the postshock shell of the matter swept up by the shock in the H-envelope. The formation of the dense shell or "wall" ) has important consequences for the evolution of the metal-rich ejecta as well as for the growth of Rayleigh-Taylor instabilities. After t pb ≈ 1 h, continuous deceleration of the forward shock causes the reverse shock to travel back in mass coordinate (see Figure 14). It takes, however, a total of approximately ∼9 h for the reverse shock to reach the inner boundary of our computational domain in model e8.8 ( Figure 13).

3D model z9.6
In model z9.6 the forward shock has already crossed the CO/He interface at t map (see Table 1) and is traveling at roughly 27, 000 km s −1 . Behind the shock, a dense shell has formed due to the deceleration of the fast neutrino-driven wind in a windtermination shock (Figure 14). In between the shocks, the high density peak atop the bulk of the ejecta is formed by the deceleration of the forward shock at the CO/He interface. Outside of the CO/He interface, with similarity to model e8.8, a featureless density profile in the He-core and hydrogen envelope leads to an untroubled but gradually decelerated expansion of the forward shock, which reaches the stellar surface at around t pb ≈ 1.3 d). At the time of shock breakout the forward shock propagates only at ∼2000 km s −1 because of its strong deceleration in the hydrogen envelope of the star. Different from the ECSN progenitor, the wind-termination shock moves inward at t pb ≈ 10-15 s (see Figure 14) to get reflected at the center before t pb ≈ 40 s. A second reverse shock forms within the He-core of the star and propagates back in radius after ∼400 s to reach the inner boundary at t pb ≈ 800 s (Figure 13), which is more than 8 hours earlier than in model e8.8.

3D model s9.0
For the 3D simulation of model s9.0 the trajectories of the forward and reverse shocks are shown in the right panels of Figure 13. The long-time simulation is initiated at t pb = 3.14 s, when the forward shock has just crossed the CO/He interface and is traveling at v sh ≈ 11, 000 km s −1 , only a fraction of the shock velocity found in the ECSN-like models. Acceleration and deceleration of the SN shock at the CO/He interface cause the formation of a dense shell in the post-shock region (see Figure 14). The density contrast between the shell and the ejecta is, however, around one order of magnitude smaller than found for the dense shells that formed at the core/envelope boundary in model e8.8, and at the CO/He interface in model z9.6. In the following, the shock slows down to v sh ≈ 7, 500 km s −1 within the He-core of the star, before it accelerates again around the He/H composition interface, reaching v sh ≈11, 500 km s −1 . This is due to the steep density gradient just below the He/H interface in the progenitor. Thereafter, the forward shock encounters the increasing ρr 3 in the hydrogen envelope and is thereby strongly decelerated, causing a compression of the postshock material into a double-peaked dense shell. At the bottom of the hydrogen shell, a strong reverse shock forms shortly afterwards (see Figures 13 and 14). Note that the formation of this reverse shock from the He/H interface occurs much later than witnessed in the ECSN-like progenitors. Eventually the reverse shock reaches the inner boundary of our numerical grid at t pb ≈ 16 h. A first reverse shock from the CO/He interface had formed around t pb ≈ 30 s ( Figure 14), but was swept outward with the expanding ejecta. The Angle-averaged shock radius R sh (black), angle-averaged reverse shock radius R rsh (red) and maximum radius of the iso-surface of mass fraction X56 Ni+Tr = 0.03 (blue). Lower panels: Average shock velocity v sh (black) and maximum velocity of the iso-surface of X56 Ni+Tr = 0.03 (blue) in our 3D long-time simulations. The horizontal or vertical gray dashed and solid lines mark the location of the He/H interface and the surface of the progenitor star, respectively. The expansion of the forward shock in the ECSN-like progenitors proceeds with continuous deceleration nearly until the shock reaches the stellar surface. In contrast, the forward shock in model s9.0 accelerates strongly around the He/H interface. The reverse shocks in models e8.8 and z9.6 form within ≈1 s and ≈100 s, respectively, whereas we witness the formation of a reverse shock from the He/H interface in model s9.0 at ≈2500 s. Note that in model s9.0 the velocity of the fastest clump containing 56 Ni+Tr is higher than the average shock velocity within the hydrogen envelope of the star. main, spherically shaped SN shock encounters the stellar surface at t pb ≈ 2.8 d, whereas the maximum radius of the shock, pushed by a giant, nickel-rich bubble, reaches the surface already at t pb ≈ 2.1 d (see discussion in Section 5.3 and Figures 21 and 22).

Morphology of neutrino-heated ejecta
In the following we focus on the long-time development of the earlytime asymmetries which we trace by the propagation of the neutrinoheated ejecta or more specifically the 56 Ni+Tr-rich material. Kifonidis et al. (2003) already noted that 56 Ni is produced between the high-entropy bubbles that expand due to strong neutrino-heating from below. Thus the distribution of 56 Ni traces the asymmetries that developed during the onset of the explosion.

3D model e8.8
In Figure 15 we show slices of the 56 Ni+Tr mass fraction of the 3D simulation of model e8.8. The dashed white line indicates the shock radius, whereas the thin dotted black line indicates the position where the enclosed mass equals the mass interior to the He/H interface in the progenitor. As in the 1D simulation, the CO and He layers are compressed into a dense shell just behind the SN shock at a radius of ≈5×10 5 km at t pb = 9 s (see Figure 14 and Figure 16). Already at t pb ≈ 1.5 s, a reverse shock begins to form at the bottom of this dense shell (see the sharp inner edge of the light-blue, narrow ring at 5×10 5 km in Figure 16 at 9 s) . By this time the bulk of the neutrino-heated ejecta has expanded in the volume inside ≈3×10 5 km. While this innermost metal-rich material seems to expand in a basically selfsimilar fashion (see similarity of the snapshots until t pb = 150 s in Figure 15), the growth of the RT instability in small protrusions at the unstable contact interface near the outer edge of the dense shell begins to corrugate the surface of the ( 56 Ni+Tr)-rich ejecta. At t pb = 300 s, the neutrino-heated ejecta of 56 Ni +Tr catch up with the expanding dense shell and are strongly decelerated. RT plumes become clearly visible near the He/H interface first at about 1400 s in Figure 15.
As a consequence, the higher entropy/lower density bubbles are compressed to flat structures, whereas the higher-density regions in between the bubbles are able to penetrate the dense shell, thereby inducing perturbations at the He/H interface. These RT plumes grow over the following hours around the mass shell of the He/H interface, mixing clumps of 56 Ni from the central volume outward into the carbon and helium rich layers.
The reverse shock, which is visible at t pb = 2.3 h in Figure 16 as the yellow to blue discontinuity, begins to propagate back in radius, thereby compressing the innermost ejecta. At t pb ≈ 10 h the inward passage of the reverse shock and the growth of the RT instability has erased the initial structures present at the onset of the explosion. However, the overall morphology of the neutrino-heated ejecta is still basically spherical with only small-scale asymmetries (see also Figure 17).

3D model z9.6
In Figure 18, we show slices of the 56 Ni+Tr mass fraction in the 3D simulation of model z9.6. The dashed white line marks the shock radius, while the dotted lines indicate the positions where the enclosed mass equals the mass interior to the CO/He and He/H interfaces of the progenitor. The dense shell that forms after the forward shock has crossed the CO/He interface (see Figure 14) is visible as the circular dark-red region at ∼1.5×10 5 km after 9 s in Figure 19. At that time the termination shock of the neutrino-driven wind can be seen at about 5×10 4 km, moving inward.
Within the first ≈9 s the fastest of the neutrino-heated ejecta encounter this dense CO-rich shell (see Figure 19) and are compressed and squeezed to flat structures around 30 s later. Similar to the results presented for the e8.8 model, the slightly over-dense regions between the high-entropy plumes induce long-wavelength perturbations as they deform and try to penetrate the RT unstable dense shell. Over the next few minutes, RT fingers grow on top of these deformations and fragment the dense shell into numerous 56 Ni-rich shrapnels.
While the fingers grow progressively, the second reverse shock from the shock propagation through the He layer (visible as the green to yellow discontinuity at t pb = 187 s in Figure 19) begins to propagate back in radius. At t pb = 829 s the reverse shock has almost reached the center of our numerical grid, having compressed and decelerated the innermost ejecta material. At t pb =1.7 h, the forward shock has crossed the He/H interface and the inner material has been fully shredded by the instability.
As the shock velocity does not change significantly at the He/H interface, we observe no additional growth of the RT instability nor the formation of another reverse shock. Thus, the morphology of the innermost ejecta seems to be determined early on, already before the shock crosses the He/H interface.
Comparing the distribution of 56 Ni and Tr of model e8.8 (Figure 15; t pb = 4.7 d) and model z9.6 ( Figure 18; t pb = 1.3 d) shortly before shock breakout, we find a slightly more clumped morphology in the 9.6 M progenitor. The structure of the neutrino-heated ejecta in model z9.6 also remains fairly spherical with many small-scale clumps (see Figures 19 and 20). However, different from model e8.8 one can recognize a hemispheric asymmetry with bigger plumes between the 8 o'clock and 10 o'clock positions in Figure 18 and weaker plumes in the opposite hemisphere. These aspherical structures go back to asymmetries that existed already in the first seconds of the explosion, visible by larger RT mushrooms and a slightly stronger shock expansion in the left hemisphere at t pb = 2 s. The strongest plume sticks out near the 11 o'clock position in Figure 20.

3D model s9.0
The evolution of the neutrino-heated ejecta in model s9.0 proceeds drastically differently from the ECSN-like progenitors (Figures 21-23). As discussed in Section 4, the initial asymmetries and shock deformation seen in model s9.0 are considerably larger than found in the ECSN-like models. Strong convection leads to the formation of a large high-entropy plume, which is rich in iron-group material and expands about two times faster than the surrounding material at the time of shock revival (t pb ≈ 0.5 s). It crosses the CO/He interface of the star at t pb ≈ 1.3 s, shortly after the forward shock. In comparison, the slowest moving material reaches the interface around 0.65 s later.
The crossing of the CO/He interface by the shock wave has several dynamical consequences. Due to the increasing ρr 3 outside of the interface, the shock is decelerated and the postshock matter is swept up and compressed into a dense shell. Note that the neutrino wind in this model is very weak and, different from model z9.6, there is no low-density central region and no wind termination shock. Instead, the central volume around the NS contains relics of ( 56 Ni +Tr)-rich low-density plumes and 56 Ni-poor, higher density downflows during the entire evolution.
The high-density shell behind the shock is aspherical in con-  Figure 14). At ≈1400 s the growing RT instability at the He/H interface begins to affect the outer layers of the neutrino-heated ejecta. From roughly 2 h on the plumes grow in size, while the reverse shock (visible at the base of the plumes) begins to propagate back in radius (see Figure 13). About 9 h after bounce the reverse shock reaches the center, gets reflected there, and on the way compresses the metal-rich ejecta in the central region. Note that the plumes are almost evenly distributed in angular direction and radial extent. Cyan colored regions in the lower right panel represent the surrounding medium embedding the progenitor. The dense shell at which the reverse shock will form (see also Figure 15) is first visible as the light yellow ring with a radius of about ≈3×10 5 km at t pb = 9 s. The panel at t pb = 2.3 h shows the formation of small RT fingers. These fingers grow with time and partly take up the neutrino-heated material (see panels at t pb ≈ 2.3 h and t pb ≈ 10 h). At t pb ≈ 4.7 d the innermost ejecta are characterized by an overall spherical shape, superimposed with the relics of the RT fingers. The reverse shock is very prominent at t pb ≈ 2.3 h and travels inward from about 3 h on (see Figure 13). It reaches the center at t pb ≈ 9 h to be reflected outward again (see lower left panel). trast to the shells found in the ECSN-like progenitors, which is also reflected by the still deformed supernova shock (see first two panels in Figure 22 for t pb = 3 and 124 s). Around the dense shell we observe the growth of large plumes, which stem from the initial asymmetries of the explosion. These can be seen in the 10 o'clock direction in Figure 21 at t pb = 124 s and in Figure 23. At the tops of these plumes small RT fingers grow, in line with the analysis of the amplification factors presented in Section 5.1. Note that at this point in time the still deformed SN shock crosses the He/H interface (see Table 1 and Figure 13). Due to the varying ρr 3 profile around this composition interface, the shock accelerates and decelerates, thereby forming a dense shell (see the yellow-green ring in Figure 22 at t pb = 742 s and the density spike in Figure 13). When the fast and dense metal-rich plume encounters the shell, it induces a high-amplitude perturbation in this RT unstable layer. From this perturbation, aided by the large initial momentum of the metal-rich plume, we observe the growth of a large RT structure, which assists the further outward expansion of the 56 Ni+Tr material of the initial, big plume. As the shock is decelerated in the H-envelope, the dense metal-rich plume retains higher velocities than the speed of the shock and thus the plume is able to deform the forward shock on its way (see Figure 22 at t pb ≥ 2.4 h).
Concurrently, the reverse shock, which forms at the bottom of the He/H interface, propagates back into the ejecta and strongly decelerates and compresses the ( 56 Ni+Tr)-rich material close to the center (see Figures 21 and 22 at t pb ≈ 7.8h), whereas the fastest, biggest plume escapes the most dramatic deceleration, although its velocity also shrinks with time (see Figure 23). The growing RT instability around the CO/He interface (clearly visible in Figure 22 at t pb ≈ 7.8h) seems to only slightly affect the outer boundary of the central ( 56 Ni+Tr)-rich material, as can be seen in Figure 21 at t pb ≥ 7.8h (in line with the small amplification factors found in this region).
Due to the strong deceleration of the innermost material, the fast plume almost fully detaches from the core material as it propagates through the hydrogen envelope of the star (see last two panels in Figures 21, 22, and 23). It encounters the surface of the star at around t pb ≈ 2.1 d, so more than half a day earlier than the spherically shaped main shock front, which reaches the stellar surface at t pb ≈ 2.8 d. Why is the large plume able to travel with such high velocities, even deforming the forward shock, while the bulk of the 56 Ni+Tr mass travels at considerably slower speed? First, the fastest 56 Ni+Tr material is, at all times, in close vicinity of the immediate postshock matter (see Figure 13). Second, after the forward shock has crossed the He/H interface, the Ni-rich plume is decelerated less than the average shock (see Figure 13), since this material is denser than its surroundings in the hydrogen layer. Third, large growth rates at the He/H interface lead to an efficient outward mixing of the dense plume within the unstable layers, and the plume can therefore also escape the strong deceleration by the reverse shock. As a consequence, the ( 56 Ni+Tr)-rich plume catches up with the forward shock in the hydrogen envelope. Due to its large momentum it deforms the outgoing forward shock in its trajectory. This is similar to a transient situation at about half an hour and 3 h in model e8.8, where the neutrino-heated ejecta in RT plumes catch up with the strongly decelerated immediate postshock material, thereby pushing the forward shock (see Figures 13 and 16).

Extent of mixing
In Figure 24, we display the normalized mass distributions of various nuclear species, including free protons from the freeze-out of NSE, hydrogen from the envelope, helium, carbon, oxygen-neonmagnesium, radioactive nickel, and the tracer nucleus of neutronrich species for all of our models at the time of shock breakout as functions of radial velocity and mass coordinate (the latter is defined as enclosed mass M(r) at radius r).
The distributions of the iron-group ejecta in the ECSN-like models (e8.8 and z9.6) have similar shapes in velocity and mass space (apart from the differences that result from the different initial progenitor composition). They are characterized by a maximum centred around 0.3×10 3 km s −1 and have a high-velocity tail which extends to ∼0.5×10 3 km s −1 (using ∆M i /M i = 8×10 −4 as a thresh- where the enclosed mass equals the mass interior to the CO/He and He/H composition interfaces of the progenitor. The initial asymmetries that develope during the onset of the explosion are still visible at t pb = 1.5 s, but they are soon compressed to flat structures as they collide with the RT unstable dense shell behind the CO/He interface (see times from t pb = 9 s to t pb = 39 s and also Figure 19). Over the next, roughly, one hour, the growing RT instability mixes the outer 56 Ni-rich layers outwards in mass coordinate in numerous small fingers. 1.7 h after bounce the RT instability has basically saturated and the first iron-rich clumps reach the He/H interface of the progenitor, where in model z9.6 no secondary RT instability occurs. The final morphology of the 56 Ni-rich ejecta has mostly lost any resemblance with the state at t map and is dominated by small-scale asymmetries. However, the overall distribution of 56 Ni and Tr remains roughly spherical with many small-scale features and only a slight global deformation, which exhibits larger and stronger iron-rich plumes between the 8 o'clock and 10 o'clock directions and weaker structures in the opposite directions. A corresponding hemispheric asymmetry is already visible at the beginning of the explosion, see top left panel for t pb = 2 s. A termination shock of the neutrino-driven wind is visible as the yellow/orange discontinuity at t pb ≈ 2 s. The shock deceleration in the He-layer leads to the formation of a dense shell that gets extremely compressed by the outward shock from the reflection of the wind termination shock at the center (see Figure 14). Over the next ≈100 s, RT plumes start to grow within the unstable layer between the two shocks, thereby fragmenting the dense shell. A second reverse shock forms when the SN shock propagates through the He-layer of the progenitor. As this reverse shock travels back in radius, similar to the results presented for model e8.8, the plumes grow to their maximal radial extent (see panels at t pb = 187 s-829 s). The final morphology of the ejecta at t pb ≈ 1.3 d resembles the late-time morphology in model e8.8. Thus, the mixing is more efficient in model z9.6 in comparison to model e8.8, although the explosion energy and the integrated growth factors are larger in the latter. Why is this the case? The answer can be found by inspecting Figures 15, 16 and 18, 19. Model z9.6 shows a larger asymmetry already at the onset of the explosion. Additionally, the density profile of the progenitor of model z9.6 exhibits less extreme declines outside of the Si/CO and CO/He interfaces than the sharp drop of the density at the edge of the degenerate core in the ECSN progenitor. This permits a shock expansion that is not quite as rapid as in model e8.8 (Figures 4 and 13). The initial asymmetry triggers more and faster growth of the RT instability at the CO/He interface in model z9.6 (see Figures 18, 19) on a time scale much shorter than the growth of the RT mushrooms in model e8.8. While in model z9.6 large RT plumes are visible already at about 100 seconds after bounce, it takes a few hours for such structures to develop at the He/H interface of model e8.8 (see Figure 15). From this result we conclude that the extent of mixing during the SN blast is not only determined by the linear growth factors, but depends strongly on the initial explosion asymmetry seeding the growth of the RT instability at the unstable composition interfaces.
In contrast to the ECSN-like models, heavy elements are mixed to large mass coordinates and velocities in model s9.0. This is facilitated by the extreme initial asymmetries at the onset of the explosion. Fast metal-rich plumes arrive quickly at the He/H interface and trigger the RT instability there, thereby transporting a significant amount of the total 56 Ni+Tr mass to large mass coordinates into the H-envelope. Outward mixing of intermediate-mass elements is driven by the growth of the RT instability which causes a fragmentation of the dense shells that form after the forward shock crosses the CO/He and He/H interfaces. Most importantly, the biggest plume is able to transport 56 Ni+Tr-rich matter to large velocities and mass coordinates. It contains about 25% of the total 56 Ni+Tr mass at t map and carries about half of that to radii well ahead of the average radius of the shock when it reaches the stellar surface. We find that about 4% of the 56 Ni+Tr-rich material travels with more than 1, 000 km s −1 and thus far ahead of the bulk of the metal-rich matter.

Long-time evolution of compact remnant properties
While the supernova shock wave travels through the progenitor star, some material falls back onto the newly formed NS (e.g. , Chevalier 1989;Woosley 1989;Zhang et al. 2008;Fryer 2009;Wong et al. 2014). This fallback is needed to explain the observed broad range of compact remnant masses (Zhang et al. 2008;Wong et al. 2014;Ertl et al. 2020). A first episode of fallback occurs when the neutrinodriven wind abates and the wind termination shock moves back towards the PNS (Arcones et al. 2007). This period happens at early times, roughly within ∼10 s after bounce. At later times fallback is driven by the reverse shocks that originate from the acceleration and deceleration phases of the SN shock passing the composition shell interfaces. As described in Sections 5.2 and 5.3, these reverse shocks propagate backward into the central volume.
In Figure 25 we show, starting from t map , the corresponding time-dependent mass accretion rate through the inner boundary of our computational grid, | M(t)|(R ib ) (top panel), the associated timeintegrated accreted mass, , and the evolution of the angular momentum of the NS, J NS (t), driven by the angular momentum that is carried by fallback material into the central volume and that is assumed to be accreted into the compact remnant (bottom panel) for all of our 3D long-time simulations. For comparison, Table 6 lists the values of the NS masses (baryonic and gravitational), total kick velocities, v tot NS , total angular momenta, J NS , and spin periods, P NS , as well as the angle between NS spin and kick vectors, θ vJ , at t map and at the end of the 3D simulations, t fin .
During the first fallback episode, i.e., during the first tens of seconds, all three models behave differently because of their dif- The seed for the biggest later structure is set already at around 3 s after bounce when the shock passes the CO/He interface. (The kink of the white dashed line for the shock surface at this time is a numerical artifact due to the shock-detection algorithm.) We first observe the growth of the RT instability in the direction perturbed by the largest and fastest initial convective plume (see, e.g., t pb = 124 s, at the 10 o'clock position), and there is less deceleration of this dense clump of 56 Ni+Tr when it penetrates the CO/He interface. It therefore begins to move far ahead of the slower 56 Ni+Tr. Shortly afterwards the plume arrives at the unstable He/H interface, inducing a large-scale perturbation there. While the slower clumps of 56 Ni+Tr are further decelerated by the reverse shock that forms at the He/H interface, the biggest clump begins to push the shock, thereby transporting a significant amount of neutrino-heated ejecta to velocities larger than the average shock velocity (see also Figure 13).

Figure 23
. 3D renderings of the X56 Ni+Tr = 0.03 iso-surface of model s9.0 at the indicated times. The color-coding represents the radial velocity of the material, and the tripod in the left panel defines the orientation of the global coordinate system. During the first roughly 800 s, the largest initial asymmetries grow to extended metal-rich plumes. Slower iron-group matter in the interior is decelerated by the collision with the dense shell formed around the CO/He interface. In the following, these fastest clumps encounter the dense shell behind the He/H and induce large-amplitude perturbations in the RT unstable layer. Consequently, large metal rich RT fingers begin to grow from the interface. At t pb = 2.8 d we find one very big and two smaller metal-rich plumes, of which the largest one penetrates the surface of the star even ahead of the average shock radius (see also Figures 21 and 22).  Table 4). The right part of the table (columns 10-16) gives the corresponding final NS properties at the end of our long-time simulations, t fin . v tot,fin NS is the final total kick velocity of the NS, J fin NS its final angular momentum, θ fin v J the angle between spin and (total) kick vectors, M fb the total fallback mass, M fin the final baryonic mass, M fin g the gravitational mass and P fin NS the spin period, adopting a NS radius of 12 km. It is assumed that all fallback matter is accreted by the NS. ferent evolution during the neutrino-wind phase. In model e8.8 the forward shock is very fast and the slower neutrino-driven wind (adopted from a parametric 1D simulation; see Fig. 3) expands freely behind the shock without deceleration and without developing a reverse shock. After the termination of the neutrino-driven wind at ∼2.5 s post bounce, fallback sets in, but the mass accre-tion onto the NS declines steeply, because the ejecta move rapidly outward, evacuating the surroundings of the compact object. In model z9.6 a similar situation applies, but the faster and more long-lasting neutrino-driven wind (adopted from a self-consistently computed NS cooling model; see Fig. 3) is stronger and evacuates the neighbourhood of the NS even more extremely than in model e8.8. However, a reverse shock forms when the forward shock passes the CO/He interface and decelerates, while at the same time the neutrino-driven wind pushes from behind and compresses the postshock matter into a dense shell. This reverse shock moves inward within a few seconds and creates the short accretion spike at ∼8 s ( Figure 25). The accretion peak decays within only a second when a reflected wave sweeps through the medium surrounding the PNS outward again. Model z9.6 displays the lowest mass accretion of all three models. In contrast, model s9.0 has by far the highest mass accretion through the inner grid boundary since its neutrino-driven wind is very weak. Therefore the NS is not surrounded by a large low-density wind bubble, but instead accretion downflows and rising plumes of neutrino-heated matter continue to coexist in the vicinity the NS for many seconds of postbounce evolution. The mass accretion rate during the early fallback phase exhibits a correspondingly high plateau between ∼3 s and ∼10 s. Models e8.8 and z9.6 reach the asymptotic scaling M∝ t −5/3 (Chevalier & Klein 1978;Zhang et al. 2008;Dexter & Kasen 2013;Wong et al. 2014) already after about 10 s, and this persists until the late fallback associated with the reverse shocks sets in at about 1000 s in z9.6 and at about 9 h in e8.8. In contrast, model s9.0 displays clear deviations from the −5/3 power-law until roughly 1 h, and only gradually approaches the −5/3 power-law behavior later, because of large-scale asymmetries in the fallback material. Late accretion due to the reverse shock from the He/H interface is triggered only after ∼16 h.
With an integral value around 5 × 10 −3 M model s9.0 has the highest total fallback mass (middle panel of Fig. 25), the other two models possess at least 10 times lower values. In all of the three cases the fallback mass is too low to have any significant impact on the NS mass or kick (see Table 6). In principle, the asymmetric fallback of matter can change the net momentum carried by the ejecta with a corresponding increase of the momentum of the NS in the opposite direction. In the considered models this modifies the NS velocity by at most a few kilometers per second (certainly 10 km s −1 ), with the biggest effect in model z9.6 because of the early onset of the reverse-shock induced accretion in this case.
However, the asymmetric fallback of matter also carries angular momentum through the inner grid boundary. Although very little mass is accreted, matter that falls back from large radii can transport appreciable amounts of angular momentum, accounting for the dominant contribution to the total NS angular momentum in all of our simulations. In model s9.0 the angular momentum of the After some initial transition phase, which lasts longer in model s9.0 because of the 3D asymmetries of the slow ejecta in the vicinity of the NS, all models adjust to the well-known power-law decline according to M∝t −5/3 (black dashed line). This continues until the reverse shock from the outer layers in each model has propagated inward and reaches the inner boundary of the computational grid at radius R ib . In models e8.8 and s9.0 early fallback clearly dominates, whereas in model z9.6 the reverse shock contributes significantly to the total fallback mass. Angular momentum carried through R ib by the infalling matter changes the initial angular momentum of the NS under the assumption that all fallback matter is accreted by the NS. In models e8.8 and z9.6 late fallback associated with the reverse shock causes the main effect on J NS , whereas in model s9.0 the total angular momentum rises continuously and most steeply at early times, when some of the highly asymmetric ejecta fall back onto the NS.
NS increases by more than a factor of 10 to 10 47 g cm 2 s −1 within the first 10 s of fallback accretion. It further grows continuously by another factor 2.5 until the 3D simulation was stopped over 4 d later. The final angular momentum of the NS is more than 30 times bigger than at 3 s after bounce. In models e8.8 and z9.6 the NS angular momentum begins to change only when the late reverse shocks reach the center and anisotopic structures begin to be accreted. This amplifies the initial angular momentum of the NS in z9.6 still by a factor of 15 (see bottom panel in Figure 25). In contrast, in model e8.8 the total effect is much more modest because of the low mass associated with the late fallback and the weak asymmetry of the explosion in this case.
Nevertheless, also in model e8.8 the spin period estimated for a NS with 12 km radius decreases from ∼11 s at 0.5 s after bounce to a final value of nearly 4 s ( Table 6). For z9.6 the spin period shrinks from 3 s at t map = 1.44 s to finally 0.19 s, and for s9.0 the final period is as low as 0.030 s, whereas it had been ∼1 s before the fallback.
The NS spins seem to be randomly oriented relative to the NS kick directions, i.e., the angles between NS spin vector and total kick vector, θ vJ , do not show any preference for spin-kick alignment, neither at early times nor days later when the fallback is complete and the simulations are terminated. This is not unexpected and it is in line with previous findings (Wongwathanarat et al. 2013;Müller et al. 2019;Chan et al. 2020). To date there is no suggestion based on well accepted physics (i.e., without invoking uncertain ingredients or extreme physical assumptions) for a convincing mechanism that could provide spin-kick alignment. Even if such an alignment were achieved during the first seconds of the explosion, when ejecta and NS are still in contact through hydrodynamical and gravitational forces and the NS is kicked by anisotropic neutrino radiation, it is very hard to imagine how this initial alignment could not be overruled by the stochastic effects of the later fallback and its dominant influence on the NS spin. It might require very rapid progenitor rotation, possibly very strong magnetic fields, to impose a preferred direction for the explosion, correlating the NS recoil acceleration and the spin-up of the NS either by angular momentum inherited during the collapse from the rotating progenitor or later through accretion of stellar angular momentum by fallback (for a suggestion of such a mechanism, see Janka 2017).
Recently, Chan et al. (2020) presented results of an explosion simulation for a non-rotating, zero-metallicity (Pop III) 12 M progenitor, including the effects of anisotropic fallback. The overall behavior of this model is similar to our model s9.0, however much more extreme, because the model has a considerably higher fallback mass (nearly 0.2 M ), and the corresponding change of the NS kick is several 10 km s −1 . The accretion of angular momentum associated with the fallback spins the NS up from an early period around 100 ms to a final period of only a few milliseconds. Interestingly, in this model as well as in the other two cases considered by Chan et al. (2020), in which black holes are formed by fallback, the final angles between spin and kick vectors are close to 90 degrees, very similar to our result for model s9.0. Chan et al. (2020) explain such a perpendicular orientation by the fact that the directions of the accretion streams remain relatively constant in time, thereby corresponding to a single, off-center momentum impulse onto the remnant. Our models e8.8 and s9.0 comply with this pattern, but model z9.6 deviates from this behavior, showing spin-kick anti-alignment.  Table 2) at the time of shock breakout. Maximum velocities and the velocities of the bulk of 56 Ni and Tr scale roughly with E exp , implying slightly enhanced mixing for more energetic models (see also Figure 27).  Table 2) at the time of shock breakout. Models with higher E exp show slightly more efficient mixing of iron-group material.

Dependence on the explosion energy
In Figure 26 we present the normalized mass distributions of chemical elements as functions of radial velocity for the 2D simulations of model e8.8, using the explosion calibrations as listed in Table 2, at the time of shock breakout.
Comparing the distributions of model e8.8 2D 10 with the 3D simulation, the extent of mixing in model e8.8 does not seem to be significantly dependent on the chosen dimensionality. This allows us to investigate the influence of the explosion energy on the efficiency of mixing in our 2D simulations of model e8.8.
Inspecting Figure 26, we find that the bulk of 56 Ni resides at low velocities of ∼(0.2−0.4)×10 3 km s −1 , increasing with explosion energy roughly as v r ∼ E exp . For larger explosion energies the downward mixing (to smaller velocities) of lighter elements also seems to be more efficient, but the effect is small and affects only 1% of the respective masses. For lower explosion energies the amount of 56 Ni+Tr experiencing fallback (v r < 0) grows.
The mixing in velocity space corresponds to a distribution of 56 Ni and Tr in mass space to a maximum mass coordinate of about 1.42−1.60 M (see Figure 27), also increasing with explosion energy, using ∆M/M = 1×10 −4 as a threshold value for plotting the distributions. Note that the 4.49 M of the hydrogen envelope of the progenitor extended initially from 1.34 M to 5.83 M , thus mixing only affects the innermost part of the envelope.
As the iron-core progenitors are exploded self-consistently, their explosion energies are fixed within our framework. Drawing direct connections between the amount of mixing and the explosion energy can therefore not be done. Comparing the z9.6 model with the electron-capture model, which both have similar explosion energies, however, suggests that the influence is secondary. More decisive for the amount of nickel mixing are the progenitor structure and initial asymmetries right after shock revival. This view is supported by the strong mixing exhibited by model s9.0, which has a comparable explosion energy as well. The velocity and density perturbations that are present at the onset of the explosion, combined with the strong acceleration and deceleration of the forward shock in the envelope of the progenitor, yield high growth rates of RT instability over a larger range in the mass coordinate and thus dominate over effects caused by different explosion energies in a given progenitor.

COMPARISON TO PREVIOUS STUDIES
Some previous works also explored the explosion properties of CC-SNe of low-mass iron and ONeMg-core progenitors. The study by Radice et al. (2017) simulated the onset of the explosion of models e8.8 n , z9.6, and s9.0 in 2D with varied microphysics. They found for model e8.8 n an explosion energy of up to 1.8×10 50 erg, which is about a factor of two higher than obtained for the same 8.8 M star by Kitaura et al. (2006);Janka et al. (2008);Hüdepohl et al. (2010);Fischer et al. (2010), andvon Groote (2014). Consistently, their 2D simulations yielded explosion energies for models z9.6 and s9.0 that were ∼50% higher than in our simulations, namely 1.2×10 50 erg and 0.7×10 50 erg, respectively; their 3D calculations of model s9.0 yielded even 1.0×10 50 erg (Burrows et al. 2020). Müller et al. (2019) focused on the onset of the explosion in helium core progenitors from binary evolution and low-mass single stars also including model z9.6. They found a slightly higher explosion energy of model z9.6 of around 1.3×10 50 erg at the end of their simulation. Moreover, their explosion energy still seems to grow fairly steeply when they stopped their simulation.
However, Müller et al. (2019) used an approximative neutrino transport treatment ("fast multi-group transport", FMT) with simplified neutrino interactions, and their explosion energies should therefore be taken with a grain of salt. The reason why Radice et al. (2017) and Burrows et al. (2020) obtained consistently more energetic explosions than in our V -P simulations is unclear to us, in particular because their FORNAX code is claimed to possess neutrino physics that is compatible with that of V -P . We note that, due to the higher explosion energies found in the studies of other groups, the final PNS masses there are smaller than ours in general. A more energetic explosion drives a stronger wind from the PNS, which carries away mass from its surface or reduces further accretion. For example Müller et al. (2019) found the PNS baryonic mass of model z9.6 to be 1.35 M , close to ours, and Burrows et al. (2019) determined the PNS mass of model s9.0 to be 1.342 M , whereas we get 1.351 M before fallback accretion and 1.356 M afterwards. The explosion energies and PNS masses of Radice et al. (2017) and Burrows et al. (2020) are not only in conflict with ours but also with other previous studies as for example Kitaura et al. (2006); Janka et al. (2008);Hüdepohl et al. (2010); von Groote (2014) for model e8.8 and Glas et al. (2019a) for model s9.0. These studies consistently attain lower explosion energies with around 10 50 erg for model e8.8 and around 0.5×10 50 erg for model s9.0, respectively.
Concerning the long-time evolution of the explosion of lowmass CCSNe progenitors, Müller et al. (2018) followed the expansion of the SN shock from its initiation by neutrino heating until shock breakout in an ultra-stripped helium star (he2.8) with a helium-core mass of 1.49 M , which is structurally similar to our model z9.6. The SN runs (s2.8) of Müller et al. (2018) and Müller et al. (2019) also explode with energies only slightly higher than our simulation of model z9.6, and comparing the maximum mass coordinate of the neutrino-heated ejecta, they exhibit a similar extent of mixing. In their simulation a small fraction of the total iron-group material is mixed to the edge of the helium core of the progenitor star, quite analogously to the case of our model z9.6, where we find a small fraction of the total neutrino-heated ejecta to be mixed out to M(r) ∼ 1.8 M .
Other studies such as those of Kifonidis et al. (2006), Hammer et al. (2010, Wongwathanarat et al. (2015), Chan et al. (2018), Chan et al. (2020), andOno et al. (2020) and Orlando et al. (2019) also performed simulations of the long-time evolution of CCSNe after shock revival, based on 2D/3D initial data of the beginning explosion. These studies, however, focused on more massive RSG progenitors or BSG progenitors that stand as a proxy for Sanduleak -69 202, the progenitor of SN1987A. Kifonidis et al. (2006) and Hammer et al. (2010) Chan et al. (2020) started with initial data from neutrino-driven explosion simulations, Ono et al. (2020) and Orlando et al. (2019) initiated the explosions by injecting energy near the IG/Si interface and parameterized the deformation of the outgoing shock wave in order to mimic the effect of non-radial instabilities at the onset of the explosion and to reproduce the observed morphology of SN1987A. Long-time 3D SN simulations for studying mixing and explosion asymmetries by Ellinger et al. (2012) and Joggerst et al. (2009);Joggerst et al. (2010) started their runs of zero metallicity, low-metallicity, and solar metallicity 15 M and 25 M progenitors from spherically symmetric explosions.
These more massive progenitors differ strongly in their ρr 3profiles when compared to our ECSN-like models. Only model s9.0 exhibits structural features (e.g., a significant variation of the density gradient at the He/H interface) similar to the RSG models of Wongwathanarat et al. (2015). Consequently, all models of the mentioned studies evolve considerably differently from the ECSN-like progenitors presented here. Larger amplification factors at the composition interfaces and a more extended region of instability lead to the growth of large RT plumes, which are absent in models e8.8 and z9.6. Model s9.0, however, behaves in a more similar way, despite the smaller ZAMS mass and considerably lower explosion energy. Although the previous studies are tuned to give around 10 51 erg for the explosion energy (e.g., to be compatible with observations of SN1978A), we find a similar efficiency of mixing in s9.0, reflected by the distributions of the chemical elements at the end of our simulation. This efficient mixing is facilitated by the strongly asymmetric onset of the explosion and the growth of strong secondary RT instability at the composition interfaces triggered by the initial ejecta asymmetries.

SUMMARY AND CONCLUSIONS
In this study we presented results of 1D, 2D, and 3D SN simulations for three non-rotating low-mass progenitors, two of which were RSGs of 9.6 M (z9.6) and 9 M (s9.0), respectively, which had formed ∼1.30 M iron cores at the end of their lives, and the third one was a newly constructed 8.8 M super-AGB star as ECSN progenitor (e8.8) with a highly degenerate ∼1.34 M ONeMg core 10 and a pre-collapse mass of 5.83 M . Our aim was a comparison of observable features between the models, including the remnant properties, ejecta composition and asymmetries, as well as the radial mixing of chemical species during the SN blast. To this end our simulations covered continuously all evolutionary phases, from the onset of stellar core collapse to core bounce, shock formation, shock stagnation, delayed shock revival by neutrino heating, shock propagation through the stellar mantle and envelope, to shock breakout from the surface of the star, and beyond this moment until fallback of matter to the central compact remnant was complete.
Our investigation by means of neutrino-hydrodynamical simulations was focussed on neutrino-driven explosions, because all of the considered progenitors explode self-consistently and fairly quickly after core bounce by the delayed neutrino-driven mechanism. We therefore did not invoke any additional effects such as "jittering jets" (e.g., Soker 2010), which have been suggested as alternative or additional mechanism to revive the stalled SN shock (e.g., Soker 2019), and whose effects have recently been claimed to play a role in the low-energy explosions of low-mass SN progenitors (Gofman & Soker 2019).

Explosion energies and NS properties
Collapse and explosion of ONeMg-core progenitors have been investigated extensively before by fully self-consistent simulations (Kitaura et al. 2006;Janka et al. 2008;Hüdepohl et al. 2010;Fischer et al. 2010;Müller et al. 2013;von Groote 2014;Radice et al. 2017), yielding quite a spread of blast-wave energies, ranging from several 10 49 erg to nearly 2×10 50 erg, depending on details such as the treatment of neutrino transport, general relativity, and the EoS of supranuclear matter in the PNS. Therefore, we simulated the explosions of model e8.8 in 1D, 2D, and 3D with the P -HOTB code by imposing suitable neutrino luminosities to tune the SN energies to values between 3×10 49 erg and 1.5×10 50 erg. Our models reproduce the generic behavior seen previously when the shock wave expands extremely rapidly at running down the steep density gradient that marks the edge of the ONeMg core.
We picked our 3D ECSN model of e8.8 to have an explosion energy of 1.0×10 50 erg, which is in the ballpark of the previous selfconsistent runs by Kitaura et al. (2006) and Hüdepohl et al. (2010), but it is only about half of the energy obtained by Radice et al. (2017). The hydrodynamic kick of the NS in this 3D model was less than 0.5 km s −1 , and in all of our 2D realizations it stayed below about 1.5 km s −1 , in agreement with the low kick values obtained by Gessner & Janka (2018). Models z9.6 and s9.0 were exploded fully self-consistently in 3D by applying the V -P code with sophisticated neutrino transport. First results of the initial ∼0.5 s after bounce for these two simulations were presented by Melson et al. (2015a) and Melson et al. (2019). The explosion energies saturate at values below 10 50 erg, in the case of z9.6 at 0.85×10 50 erg, and for s9.0 at about 0.5×10 50 erg (compatible with Glas et al. 2019a). This is again roughly a factor of two lower than the energies found in 2D and 3D simulations for s9.0 by Radice et al. (2017), Burrows et al. (2019), and Burrows et al. (2020), and in 3D for z9.6 by Müller et al. (2019). However, all of our 3D models explode with energies consistent with the 5×10 49 -10 50 erg of the Crab SN (Yang & Chevalier 2015), which has been interpreted as ECSN or as CCSN of a low-mass ironcore progenitor (e.g., Smith 2013;Tominaga et al. 2013). Moreover, all models produce small ejecta masses of radioactive 56 Ni, beween about 4×10 −3 M and roughly 6×10 −3 M .
Both iron-core SN models develop a pronounced dipole mode of the lepton-number emission by neutrinos, which goes back to the LESA (Lepton-Emission Self-sustained Asymmetry) phenomenon. The NS kicks induced by the corresponding dipole component of the total neutrino luminosity are around 25 km s −1 , which is subdominant compared to the hydrodynamical NS kick (∼30 km s −1 ) in s9.0 but higher than the hydrodynamical kick (∼10 km s −1 ) in z9.6. The neutrino-induced kicks might further increase after the end of the evolution that we simulated with detailed neutrino transport (about 0.5 s after bounce), although in model z9.6 a quadrupole emission mode begins to dominate the dipole mode at that time, and the NS acceleration by anisotropic neutrino radiation becomes correspondingly small.
All of our SN simulations of the low-mass progenitors produce NSs with baryonic masses between ∼1.30 M and ∼1.35 M , corresponding to gravitational masses between ∼1.20 M and ∼1.23 M . Right after their formation, during the first seconds of their lives, the new-born NSs spin slowly with periods in the range of seconds, corresponding to angular momenta of the order of some 10 45 g cm 2 s −1 . However, this initial spin is dwarfed by later fallback effects. Despite the small fallback masses of at most a few 10 −3 M , which neither change the NS masses nor NS kicks to any relevant extent, the fallback transports large amounts of angular momentum to the compact remnant. The angular momentum received from fallback material outruns the previous angular momentum by up to a factor of 30 and can shrink the NS spin period from seconds to tens of milliseconds.
For example, the compact remnant in model s9.0 thus attains a spin period of nearly 1 s before fallback due to asymmetric mass ejection during the early stages of the explosion that transferred an angular momentum of only 8×10 45 g cm 2 s −1 to the NS. After fallback the spin period of the new-born NS is 30 ms only. This is very close to the 17-19 ms estimated for the birth period of the Crab pulsar in SN 1054 (Manchester & Taylor 1977;Bejger & Haensel 2003;Lyne et al. 2015). The NS in this model received a kick velocity of 41 km s −1 . However, the velocity is still rising when we stopped monitoring the neutrino-(LESA-) induced component of the NS kick (reaching ∼25 km s −1 in our models), and the final value of the total kick velocity might well be higher. Therefore, we consider this result to be in the ballpark of the spatial velocity of the Crab pulsar, which is inferred to be around 160 km s −1 with rather big uncertainties (Hester 2008;Kaplan et al. 2008). A possible spin-kick alignment of the Crab pulsar (see the detailed discussion in Kaplan et al. 2008), however, cannot be explained by our model; s9.0 yields a final angle of ∼100 • between NS spin and kick. Since the dominant mechanism for spinning up the NS is fallback, it is very difficult to understand how such a late-time effect, whose angular momentum is connected to stochastic asymmetries of the fallback matter, could correlate with the NS kick direction, which is determined in the first seconds after the onset of the explosion (see also Chan et al. 2020, for similar conclusions for explosions of more massive progenitors). It is also not easy to imagine a scenario where either rotation of the progenitor star (which is not taken into account in our pre-collapse models) or the inclusion of the NS motion in the modelling (which we do not follow because of our use of an inner grid boundary for the long-time runs) could lead to spin-kick alignment as a deterministic consequence of physical effects. We therefore hypothesize that the spin-kick alignment of the Crab pulsar, if true and a relic of the SN explosion and not just a projection effect, is a purely incidental outcome.
From all of the three progenitors, only s9.0 seems to have favorable properties to explain the explosion energy of SN 1054 as well as the magnitude of the spin and kick of the Crab pulsar. Models e8.8 and z9.6 explode too symmetrically and possess too little fallback to yield the short spin period. Moreover, in such symmetric explosions the NS kick is strongly dominated by a component associated with anisotropic neutrino emission (which our simulations tracked only in models z9.6 and s9.0). Model s9.0 also demonstrates that very low explosion energies (∼0.5×10 50 erg in this case) do not exclude sizable NS kick velocities when the explosion occurs highly asymmetrically; the momentum asymmetry of the ejecta is α ej ∼10% after ∼3 s in this model. Therefore we reason that the observed kick and spin of the Crab pulsar may be most naturally explained by considerable asymmetries of the explosion and fallback in a CCSN of a low-mass Fe-core progenitor. It is also important to note that the considered progenitors of z9.6 and s9.0 are just two samples of this class of SN progenitors, and there is a large variety of them, filling the gaps in between, and beyond, including more extreme hAGB-like models beyond z9.6. Therefore our study is just a starting point in exploring this most interesting regime, which is not small in terms of its weight by the stellar initial mass function (e.g., for the z-series it spans the range of ∼9.6-10.3 M ).
Nevertheless, alternative possibilities cannot be excluded on grounds of our results, because our conclusions apply for the considered low-mass progenitor properties, which do not include rotation at the onset of core collapse. Rotation of the progenitor core can be relevant for the spin of the NS, in particular in explosions with little fallback or little angular momentum connected to the fallback (for a recent discussion of the many facets of angular momentum transport in massive-star evolution and possible implications for NS rotation, see Ma & Fuller 2019). For example, the degenerate pre-collapse core may spin up considerably during contraction, and the compact ONeMg core near the Chandrasekhar mass limit might rotate rapidly, if the initial angular momentum is (partially) conserved. This latter requirement may be enabled by the fact that the core evolves very much independently from the envelope due to the very steep density gradient around the core-envelope interface of a super-AGB star. Thus, the angular momentum transport from the core to the extended envelope could be small. Rotation could even increase the mass of the degenerate ONeMg core compared to the non-rotating case ("super-Chandra" cores), with a mass excess that depends on the degree of differential rotation (Uenishi et al. 2003;Benvenuto et al. 2015;Hachisu et al. 2012).
Another possibility could be a close binary progenitor scenario, in which, for example, the formation and collapse of an ONeMg core might occur as a result of the merging of two white dwarfs in a certain common-envelope configuration formed during the close binary evolution (Nomoto 1985). In such a case, the disrupted white dwarf material would form a relatively dense envelope around the collapsing ONeMg core and affect the properties of the resulting NS. This scenario as well as the previous one might account for the spin of the Crab pulsar, but both of them are likely to share the problem with our explosion models of e8.8 and z9.6 that the NS kicks stay too low.
If the NS progenitor was a member of a close binary system, as was speculated for the Crab pulsar (Tsygan 1975), the NS properties may have been affected by the binary nature. In this case the NS kick might originate from the breakup of the binary when SN 1054 exploded (Blaauw 1961) instead of an intrinsic SN kick. However, spin-kick alignment of the Crab pulsar in such a scenario would require extreme fine tuning of the binary evolution (Horvat et al. 2018).

Ejecta composition, asymmetries, and mixing
Our three 3D models exhibit considerable differences in their ejecta morphology and long-time evolution. The degree of asymmetry and extent of radial mixing show a clear dependence on the steepness of the density profile around the degenerate core. If the density drops steeply, as in the ECSN-like models of e8.8 and z9.6, quick shock revival and fast shock expansion favor buoyant plumes and bubbles from postshock convection to freeze in on relatively small scales when the shock starts to accelerate outwards. Consequently, asymmetries in the ejecta possess small scales, corresponding to higher-order spherical harmonics modes. In contrast, in model s9.0, which has a significantly flatter density decline around the iron core and higher mass-accretion rate during the shock stagnation phase, the ∼300 ms delay for the onset of the explosion permits the development of large-scale asymmetries with dominant dipolar and quadrupolar deformation modes.
Despite these differences, both iron-core models resemble each other in the directional Y e variations imposed on the neutrino-heated material by the neutrino-emission dipole of the LESA: neutron-rich material is predominantly ejected in one hemisphere, whereas an excess of proton-rich matter is expelled on the opposite side. The width of the mass distribution, however, is broader in the case of z9.6, with Y e reaching down to nearly 0.39 and up to ∼0.63, in contrast to 0.46 Y e 0.58 in s9.0.
In both the e8.8 and z9.6 progenitors, the SN shock expands very rapidly at the beginning, but is dramatically decelerated when it travels through the hydrogen plus helium shells with their positive gradients of ρr 3 . This triggers the formation of a strong reverse shock propagating inward as well as the creation of conditions for the growth of RT instability, which distributes the neutrino-processed material including freshly nucleosynthesized radioactive species like 56 Ni in stretched fingers and plumes within an extended spatial volume. Nevertheless, very little of this material gets mixed into the hydrogen envelope, and the corresponding distribution remains narrow in mass space, stretched out only over the innermost ∼0.1 M of the ejecta in the case of e8.8 and ∼0.6 M in z9.6. In both cases the maximum 56 Ni velocities are around 500 km s −1 .
In contrast, in model s9.0 the SN shock propagates rather steadily, yet less rapidly (though still with ∼10 4 km s −1 ), during the first 150 s, but it is also strongly decelerated after entering the hydrogen envelope. Because the shock is highly deformed in this explosion and the postshock ejecta are extremely asymmetric from the beginning, these initial ejecta asymmetries trigger the rapid and powerful growth of Richtmyer-Meshkov instability and RT instability at the He/H interface. Big plumes of nickel-rich matter, originating from the biggest bubbles at the time of shock revival, shape the large-scale asymmetry of the SN blast by penetrating deep into the hydrogen envelope. In fact, the most extended of these plumes is pushed by buoyancy forces so strongly that it catches up with the continuously decelerated SN shock and creates a massive deformation of the shock front, thus overtaking the average shock radius. This effect persists until the shock breaks out from the stellar surface, which therefore happens highly asymmetrically. The first, biggest plume pushes the shock through the surface of the progenitor after roughly 2.1 days, whereas the main sphere of the shock reaches the stellar surface considerably later after 2.8 days. Since radioactive nickel is mixed through the entire hydrogen envelope in this model, it expands with velocities up to 1400 km s −1 after the breakout from the star (and might be even further accelerated due to radioactive decay heating).
The bubble-driven, asymmetric breakout of the SN shock in our 9 M model will have ramifications for theoretical studies of this evolution phase and for the interpretation of corresponding observations (e.g., Nakar & Sari 2010). Since at that time the giant plume contains about 10-20% of the radioactive material produced in the deepest layers of the SN, the one-sided expansion of such a feature will cause a strongly direction-dependent emission of gamma-rays and X-rays, which will occur particularly early in the hemisphere of the bubble breakout (for recent investigations of high-energy radiation emission based on asymmetric 3D explosion models, see Alp et al. 2018Alp et al. , 2019Orlando et al. 2019;Jerkstrand et al. 2020). The deep mixing of iron-group nuclei and radioactive species from the core through the whole hydrogen shell also plays an important role in shaping the Type-II SN light curve during the luminosity peak (Utrobin et al. 2017). Moreover, the large-scale deformation of the stellar envelope by extended, wide-angle plumes enriched with heavy elements might also offer an explanation why some Type-IIP SNe show an unusually early rise of the polarization before the tail phase is entered, and thus before the helium core is exposed by the transparency of the hydrogen envelope (Nagao et al. 2019). Last but not least, an outward pushing RT plume with its long-stretched stem might be a mechanism to create the faint protrusion extending out from the northern rim of the visible Crab Nebula, which is often called northern ejecta 'jet' (e.g., Blandford et al. 1983;Davidson & Fesen 1985;Fesen & Staker 1993;Black & Fesen 2015). Such an origin would naturally allow for understanding the fact that the jet's sharp western limb and its blueshifted and redshifted sides seem to be radially aligned with the center of expansion of the remnant (Black & Fesen 2015). The hollow appearance of the jet in [OIII] line emission, its remarkably empty, elliptical shape, and its growing diameter with larger distance from the remnant center, all of which define a funnel-like structure, could be naturally explained in such a picture. The jet walls would be expected to contain oxygen swept up when the RT plume penetrates the oxygen shell of the progenitor, whereas the jet's interior should contain a reduced oxygen fraction but enhanced content of iron-group elements and possibly also silicon.
A more detailed analysis of the evolving plume structure, based preferably also on higher-resolution 3D simulations, is needed to consolidate such an appealing scenario. Another interesting extension of our work concerns the nucleosynthetic post-processing of the LESA-affected neutron-rich and proton-rich components of the neutrino-heated ejecta. Since our 2D and 3D results of the mass distributions as functions of Y e exhibit, overall, fairly close similarities, however, we expect basic confirmation of the trends already seen in a recent study of element formation in the ejecta of 2D simulations for low-mass stars including the e8.8 and z9.6 progenitors (Wanajo et al. 2018). Finally, it would be desirable to repeat our 3D explosion modeling with 3D pre-collapse conditions originating from the latest phases of convective oxygen and silicon shell burning (e.g., Arnett & Meakin 2011;Couch et al. 2015;Müller 2016;Müller et al. 2016;Yoshida et al. 2019;Yadav et al. 2019). The large-scale density and velocity perturbations created by the convective shell burning prior to collapse are more realistic seeds for the growth of postshock instabilities than the artificial, small-amplitude, stochastic cell-by-cell seed perturbations that we imposed on the spherical progenitor models in our present calculations. The larger physical perturbations of 3D initial models might lead to earlier explosions, in particular of the s9.0 progenitor, and might also affect the shock and explosion asymmetry right after shock revival (see Couch & Ott 2013;Müller et al. 2017;Kazeroni & Abdikamalov 2019).

APPENDIX A: PNS COOLING MODEL AND INNER BOUNDARY CONDITION IN PROMETHEUS-HOTB
As stated in Section 3.4, we use the modeling approach of Ugliano et al. (2012), Sukhbold et al. (2016), and Ertl et al. (2020) in simulations with P -HOTB. The central 1.1 M of the PNS are excised from the computational domain and replaced by an inner grid boundary at R ib . The shrinking of the cooling and deleptonizing PNS is mimicked by the contraction of the inner grid boundary, whose time dependence (see also Arcones et al. 2007) is given by where R ib,f is the final radius, R ib,i the initial radius and t 0 the contraction timescale. R ib,f and t 0 are two representatives of our set of free parameters and are chosen to mimic the behavior of the PNS contraction found in more sophisticated simulations of PNS cooling (see Scheck et al. 2006;Sukhbold et al. 2016, for comparisons with such results). As detailed in Ugliano et al. (2012), the PNS core of mass M c = 1.1 M is described by an analytic one-zone model under the constraints of energy conservation and the virial theorem including the effects associated with the growing pressure of the accretion layer. The accumulation of mass around the PNS core is followed by the hydrodynamic simulations. The one-zone model provides the time-dependent total neutrino luminosity that leaves the excised core and is imposed as boundary condition (split up into contributions of each of the neutrino species) at the bottom of the computational domain at R ib . It is given by Here, Γ = 3 is the adiabatic index of the PNS core (assumed to be homogeneous), G is the gravitational constant, M c is the core mass, a is a parameter which characterizes the accretion luminosity, and m acc (t) is the mass contained between the radius of the inner grid boundary, R ib (t), and the radius r 0 at which the density falls below ρ 0 = 10 10 g cm −3 . We define m acc (t) = − 4πr 2 0 v 0 ρ 0 , where v 0 is the fluid velocity at the position r 0 . In multi-dimensional simulations, m acc is determined from angle-averaged values. The time dependence of the core radius R c in Equation (A2) is prescribed by where p < 0 is another parameter. We always set the characteristic time scale t L = 1 s and use initially R c,i = R ib,i . Note that in general the PNS core radius and the radius of the inner grid boundary can differ during parts of the evolution. The quintuple of p, R ib,f , a, R c,f , and t 0 constitutes our set of five parameters to approximate the physics of the time evolution of the PNS and to enable supernova explosions with chosen energy. The calibration of these parameters was done using the method described in Ertl et al. (2016).

APPENDIX B: CORRECTION TO NEUTRINO-NUCLEON SCATTERING IN PROMETHEUS-HOTB
Because of numerical issues in the neutrino transport module of the previous version of P -HOTB, long-time simulations (t pb > 3 s) including our neutrino transport approximation showed spurious oscillations in the energy source terms Q ν . In particular, cases where no or only a very late explosion was observed were affected. These undesired effects were caused by an improper treatment of the energy source terms connected to our non-conservative description of neutrino-nucleon scattering. Scheck et al. (2006) coined the net energy exchange rate through neutrino-nucleon scattering, following Tubbs (1979) and Janka (1991), in a closed form: which is their equation (D.68). Here, the n (n being a natural number) are spectal averages of powers of the neutrino energy, and T is the gas temperature in MeV. When the temperature exceeds 4 /6 3 , the net rate Q νN is negative and neutrinos receive energy from the stellar medium through scattering reactions with nucleons. Since the scattering term was implemented in this closed form and could change its sign, it appeared either as an opacity producing contribution to the neutrino-energy absorption rate Q − (when Q νN was positive) or as an energy source rate for neutrinos, Q + (when Q νN was negative), in the analytic integral of the transport equation employed by Scheck et al. (2006). Because of the associated big changes and alternating signs, the tight coupling of fluxes and source terms caused large-amplitude oscillations in the transport solution for the neutrino fluxes, induced by large variations of the heating and cooling source terms Q + and Q − . This led to unphysical heating of the matter, creating additional luminosity and preventing further cooling of the PNS.
A simple solution to this problem is to split Equation (B1) into two separate source-/sink terms, a neutrino-energy emission rate represented by the temperature-dependent term, and a neutrino-energy absorption rate, The latter rate corresponds to an exponential attenuation factor of the luminosity connected with the neutrino absorption opacity. This factor also accounts for the re-absorption of neutrinos that are locally produced in each computational cell, and thus it can damp source-rate variations very efficiently.

APPENDIX C: MULTIPOLE DECOMPOSITION OF THE NEUTRINO LEPTON-NUMBER FLUX
For our discussion, we decompose the electron-neutrino leptonnumber flux density F n ν e − F n ν e into spherical harmonics. F n ν e and F n ν e denote the individual radial number flux densities of ν e andν e , respectively. The real spherical harmonics are defined as with normalization factors and associated Legendre polynomials P m (cos θ). The coefficients for the multipole analysis are In our chosen normalization, the reconstruction reads From the coefficients c m , the multipole moments can be calculated according to Tamborra et al. (2014b) used a different definition for the monopole and dipole components of the lepton-number flux density. Their A Monopole is equal to our A 0 , whereas A Dipole = 3A 1 . This can be easily seen if we consider the lepton-number flux density to consist only of monopol and dipole components. If we also align the coordinate system into the dipole direction, the reconstruction of Equation (C4) reads Note that c −1 1 and c 1 1 vanish in this orientation of the coordinate system. Expressed in terms of the multipole moments, the latter expression becomes According to Tamborra et al. (2014b), F n ν e (θ, φ) − F n ν e (θ, φ) is proportional to A Monopole + A Dipole cos θ, if higher-order multipoles do not contribute. We can directly infer A Dipole = 3A 1 .

APPENDIX D: NEUTRON STAR KICK AND SPIN
Two mechanisms are considered that can lead to a recoil kick of the newly formed NS during the SN blast. First, asphericities developing during the explosion exert hydrodynamic and gravitational forces on the PNS. Both accelerate the PNS in the direction opposite to the strongest direction of the explosion, compatible with global momentum conservation. Since the gravitational effects become dominant during the long-time evolution, this kick mechanism was termed "gravitational tug-boat mechanism" (Wongwathanarat et al. 2013). For a thorough discussion of its physics details, the reader is referred to Scheck et al. (2006), Wongwathanarat et al. (2013), Janka (2017), Gessner & Janka (2018), and Müller et al. (2019). Second, anisotropic emission of neutrinos, specifically the anisotropic neutrino energy flux density, F e ν (θ, φ), radiated from the surface of the PNS (θ and φ are the direction angles in a polar grid of the star), exerts a force onto the PNS in the direction opposite to the most intense neutrino emission.
Using the momentum conservation equation, the hydrodynamic PNS kick can be simply estimated as where v v v hyd NS is the hydrodynamic kick velocity, M b is the baryonic (PNS) mass contained inside the radius R NS where the angleaveraged density drops below 10 11 g cm −3 , and P P P gas = ∫ R sh

R gain
ρv v vdV is the total linear momentum of the ejecta between the gain radius, R gain , and the SN shock, R sh . The momentum transfer by escaping neutrinos is given by where F e ν is the neutrino energy flux summed over all species, c the speed of light, e e e r the unit vector in radial direction, and R free the radius of evaluation (typically about 400 km), exterior to which neutrinos stream essentially freely and a tiny fraction of still interacting neutrinos can be ignored. Using P ν (t) = ∫ t 0 P ν (t )dt , the total kick velocity of the PNS at any time t can be calculated as v v v tot NS (t) = P P P hyd NS (t) + P P P ν NS (t) M b (t) = − P P P gas (t) + P P P ν (t) One can characterize the asymmetry of the ejecta and neutrino emission by means of anisotropy parameters. The hydrodynamic parameter reads α ej = |P P P gas | P ej , where is the total momentum stored in the ejecta, which becomes equal to the total radial momentum when the ejecta expand essentially radially. For the neutrino anisotropy parameter we use the total energy loss rate in neutrinos, which is given by The time-dependent total flux of neutrino momentum through the sphere at R free is given by c −1 E ν allowing us to define the instantaneous neutrino emission anisotropy parameter as In analogy to the linear ejecta momentum at a time t, the momentum radiated by neutrinos until time t is so that the time-integrated neutrino emission asymmetry at time t becomes In addition, we compute the PNS spin by integrating the flux of angular momentum through a sphere of radius r 0 around the origin, where ρ is the matter density, v r the radial velocity, and r r r and v v v are the position and velocity vectors, respectively. During the early postbounce evolution with neutrino treatment, r 0 = 100 km, whereas we use r 0 = R ib during the long-time simulations. In order to estimate the spin period of the PNS P NS = 2πI NS /|J NS |, with I NS being the moment of inertia of the PNS, we use the approximation by Lattimer & Schutz (2005), where M g,M is the gravitational mass of the PNS in units of M , and R NS,km is the radius of the PNS in units of km. In the above equation, the gravitational mass M g can be estimated from the baryonic mass M b as (Lattimer & Prakash 2000) M where β = GM g /R NS c 2 .

APPENDIX E: SIMPLIFIED NEUTRINO TREATMENT IN VERTEX-PROMETHEUS
The computational demands of the full-fledged neutrino transport in our V -P code are a severe obstacle for continuing simulations well beyond post-bounce times of ∼ 0.5 s. In order to follow explosions to much later times including neutrino effects, we newly implemented a simplified neutrino treatment based on a lightbulb-like scheme for neutrino emission and absorption. Within this framework, we do not solve the neutrino transport equations (i.e., we switch off the transport module V ). Instead, we obtain local neutrino source terms for application in the hydrodynamics module P from an analytical scaling and transformation of the neutrino source terms that are present in a model at the end of the V -P simulation. This moment in time is chosen to be in an evolutionary phase where the model is well on the way to a successful explosion.
Doing so, we aim at capturing the most crucial neutrino effects after the onset of the explosion in a computationally efficient way, while keeping numerical transients at a minimum and thus ensuring a seamless continuation of our simulations after switching off the V transport. To further increase the time step, we remap the 3D (hydrodynamical) PNS data within 10 km to 1D and slightly reduce the radial resolution in the PNS interior at the beginning of our simulations with simplified neutrino treatment. Moreover, adding more radial zones, the outer grid boundary is shifted from initially 10.000 km to ∼ 70.000 km to be able to follow the outward propagation of the shock through the exploding star at later times (t 1 s after bounce). In the following, we will elaborate on our new approximate neutrino treatment.

E1 Source terms for energy and lepton number
In our new simplified approach, we apply the following expressions for the net neutrino cooling and heating rates per volume (i.e., the energy source terms): where Equation (E1) is employed in regions where Q 0 erg (x) < 0 (i.e., PNS cooling), while Equation (E2) is applied where Q 0 erg (x) > 0 (i.e., gain-layer heating). 11 The quantities with superscript or subscript "0" are the angle-averaged net heating and cooling rates, Q 0 erg , 11 We use Equation (E2) only outside of the PNS, i.e. at radii r > R NS , while we do not allow for heating inside the PNS. Analogously, we do not allow for an increase of the electron fraction by neutrino sources via Equation (E6) in the interior of the PNS.
the angle-averaged density, ρ 0 , and the angle-averaged temperature, T 0 , at time t 0 when we switch from the transport calculation with V at t t 0 to our simplified neutrino treatment at t > t 0 . 12 To adjust the radial profile of heating and cooling to the contraction of the gain radius, which roughly follows the contraction of the PNS radius, we define the variable x ≡ r 0 /R gain (t 0 ) = r(t)/R gain (t), which connects the radial coordinate r = r(t) with r 0 = r(t 0 ). The factor [R gain (t 0 )/R gain (t)] 2 results from the transformation responsible for the inward shift of the heating profile and ensures that the net heating rate drops like r −2 at large radii.
The functional ansatz of Equations (E1) and (E2) fulfills the requirement of a continuous and smooth transition of cooling and heating before and after t 0 . The scaling of the heating and cooling rates is motivated by the rough scaling of the ν e andν e absorption rates with the density of a gas of free nucleons, and the scaling of (nondegenerate) electron and positron capture rates on nucleons with the temperature (see, e.g., Janka 2001). The parameter a in Equation (E1) is chosen to be of order unity, depending on the model under consideration, and adjusted to best reproduce the PNS contraction behaviour obtained in a corresponding 1D PNS-cooling simulation with mixing-length convection. For the cases considered in this work, we found a = 2 to be a reasonable choice. The factors f L and f E in Equation (E2) contain the dependencies of the neutrino heating rate on the radiated luminosities L ν e and Lν e (for the energy transferred by ν e andν e absorption per unit of time) and on the mean squared neutrino energy (for the basic energy dependence of the absorption cross section), which we replace by the squared arithmetic average of the mean energies E ν e and Eν e of ν e andν e leaving the PNS. In the spirit of our approach explained above, we therefore use the following functional scaling prescriptions: L ν e (t) + Lν e (t) L ν e (t 0 ) + Lν e (t 0 ) , which holds for core-luminosity-dominated conditions, and Guided by the findings of Müller & Janka (2014), f E is assumed to scale with the square of the (baryonic) PNS mass, M b . The approximate time dependence of the ν e andν e luminosities and of the corresponding mean energies is adopted from the PNS-cooling behavior in a 1D explosion simulation of the z9.6 progenitor with V -P , which results in a NS of similar mass (M b = 1.36 M , M g = 1.26 M ) as obtained in our 3D simulations of the z9.6 and s9.0 models (see Table 4). The mass scaling in Equation (E4) shall account for a possible difference in the evolution of the PNS mass between the 1D and 3D models.
In a manner analogous to Equations (E1) and (E2), we apply the following source terms for the net rates of change (per volume) of the electron number density by neutrino emission and absorption: Q + lep (r) = Q 0 lep (x) · ρ(r) ρ 0 (r 0 ) 12 The radial profiles of Q 0 erg , ρ 0 , T 0 , and Y 0 e of Equation (E5) are smoothed by time-averaging over the last few ms of the calculations with detailed transport.
where Q 0 lep denotes the angle-averaged net change-rate of the electron density and Y 0 e the angle-averaged electron fraction at time t 0 , when we switch from the calculation with detailed neutrino transport to our simplified neutrino treatment. Deleptonization, as per Equation (E5), applies in regions where Q 0 lep (x) < 0 (and Y e > 0.01), whereas for Q 0 lep (x) > 0 (and outside of the PNS, cf. footnote 11), the electron number density is changed according to Equation (E6). Since outflowing ejecta in the neutrino-heating layer gain electron number by neutrino absorption while accretion flows deleptonize, which our angle-averaged source terms cannot capture in detail in multi-dimensional conditions, we include a factor f Y e in Equation (E6). It is defined in terms of the electron fraction (Y e ) as with v r denoting the radial velocity and Θ being the Heaviside step function. The negative sign of Equation (E7b) is supposed to approximately account for the fact that accretion downflows settling onto the PNS lose electron number. Applying Equation (E7a) shall ensure that Y e in the ejecta is limited by its kinetic equilibrium value of Y eq e = 1 + λν e /λ ν e −1 , which is approached when ν e absorption on neutrons is balanced byν e absorption on protons (Qian & Woosley 1996). For the capture rates of ν e andν e , λ ν e and λν e , respectively, we employ equations (5)-(8) of Pllumbi et al. (2015), which provide the detailed expressions including corrections associated with the neutron-proton mass difference and with weak magnetism (Horowitz 2002). The time-dependent neutrino luminosities and energy moments needed to evaluate the neutrino absorption rates are again taken from a corresponding 1D simulation of PNS cooling. 13 The gain radius in Equations (E1), (E2), (E5) and (E6), which governs the radial migration of the heating and cooling profile, also needs to be prescribed as a function of time, because its behavior can be tracked only with self-consistent neutrino transport. In our simple heating and cooling treatment we couple the evolution of R gain (t) to the contraction of the PNS radius R NS (t) in the following way: where C 0 = R gain (t 0 )/R NS (t 0 ) denotes the ratio of the gain radius to the PNS radius at time t 0 , and b = min{1.01, C 0 } determines the assumed asymptotic value of R gain (t) at late times when the mass accretion rate M(t) has declined to an insignificant level. This prescription accounts for an inflated PNS mantle and therefore increased gain radius due to ongoing accretion, and it ties in continuously with the evolution of R gain (t) at times before we switch our neutrino treatment. M is evaluated at a fixed radius of 100 km and for downflows (v r < 0) only. A note of caution is indicated here. We noticed that our cooling prescription according to Equation (E1) can lead to runaway cooling, because a built-in mechanism of self-regulation is missing 13 Constraining Y e from the high side in Equation (E7a) is relevant only when our approach is applied in spherically symmetric explosion models. In our 3D simulations we found that Y e always remains fairly close to 0.5 or the equilibrium value because of the effects of coexisting outflows and downflows. Artificial regulation towards the equilibrium value by the factor f Ye in Equation (E7a) did therefore not become active in the 3D cases.
in the cooling rate. This can affect regions in the deep interior of the PNS at late evolution times (t 1 s after bounce), but it can also cause artifacts already early on in progenitors that are more massive than the ones considered in the present work. In massive stars it can happen that powerful accretion flows reach down to, or even below, the gain radius, in course of which cooling in the narrow layer between the neutrinosphere and the gain radius is overestimated, causing accelerated PNS contraction. An optical-depth-dependent exponential damping term turns out to be ineffective in such a situation, because the artificial effect occurs in a region of rather low optical depth. Devising a cure of these problems is currently work in progress. In the simulations of z9.6 and s9.0, the cooling implementation described above worked well and ensured a smooth and continuous contraction of the PNS with a gradient nicely extrapolating the behavior of the V -P simulations. Nevertheless, in the case of the s9.0 progenitor we encountered the mentioned issues with runaway cooling in a thin shell in the deep interior of the PNS. Instead of fixing them by a hard cut-off of the local cooling, we globally switched off the anyway low energy and lepton-number loss rates of Equations (E1) and (E5), respectively, at ∼ 1 s after bounce and continued our simulation with gain-layer heating only. Since the PNS radius at that time had decreased to 20 km already and its further contraction was only slow, we could not witness any transient or discontinuous behavior as a consequence of our measure. In order to account for the subsequent shallow PNS contraction seen in the guiding 1D PNS-cooling simulation, which would shift R gain in the 3D run closer in (see Equation (E8)) and thus would enhance the heating of the PNS surroundings, we scale Equations (E2) and (E6) with another factor Here, R 1D NS is adopted from the corresponding 1D model that also provides the time-dependent neutrino luminosities and mean energies used in the scaling factors of f L and f E in Equations (E3) and (E4), as well as in Y eq e . This scaling increases the source rates Q + erg and Q + lep when the PNS radius in the 1D simulation with V -P becomes smaller than in our 3D simulation with simplified neutrino treatment, thus ensuring that neutrino heating at the bottom of the SN ejecta is not underestimated.

E2 Neutrino pressure correction
Switching off the neutrino transport would also lead to a sudden drop of pressure support in the high-density regime, if the contributions of neutrinos to the total pressure were ignored. To avoid unphysical artifacts (such as PNS oscillations), we replace the neutrino-momentum source term in the hydrodynamics equations by an adequate neutrino-pressure contribution that is added to the gas pressure.
Assuming that neutrinos are in local chemical equilibrium with matter at sufficiently high densities, we can employ the following analytic expression (using the sums of relativistic fermi integrals from Bludman & van Riper 1978): p eq ν = 4π(k b T) 4 3(hc) 3 21π 4 60 + 1 2 η 2 ν e π 2 + where T is the local temperature, k b Boltzmann's constant, and η ν e = µ ν e /(k b T) the degeneracy parameter of electron neutrinos,