EDP Sciences
Free Access
Issue
A&A
Volume 552, April 2013
Article Number A126
Number of page(s) 25
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/201220636
Published online 15 April 2013

© ESO, 2013

1. Introduction

Young pulsars have been measured to possess high space velocities with an average value of several 100 km s-1 (e.g., Arzoumanian et al. 2002; Hobbs et al. 2005; Faucher-Giguère & Kaspi 2006). Some of them move with an estimated transverse velocity even higher than 500 km s-1, such as the compact stellar remnant RX J0822-4300 in Puppis A with a velocity of 672 ± 115 km s-1 (Becker et al. 2012), and PSR B1508+55 with a velocity of  km s-1 (Chatterjee et al. 2005). Such velocities are much too high to be explained by the disruption of close binary systems in the supernova (SN) event that gave birth to the neutron star (NS). Further evidence of natal kicks and associated NS spin-up can be deduced from orbital and spin parameters and kinematic information of double NS systems (e.g., Wong et al. 2010; Farr et al. 2011; Lai 2001).

Numerous models for producing the kicks have been proposed, either attempting to explain the NS recoil by anisotropic neutrino emission during the Kelvin-Helmholtz cooling of the nascent remnant or by asymmetric mass ejection in the SN explosion (cf. Lai et al. 2001, for a review). Both scenarios constitute fundamental possibilities that exploit momentum conservation in the transition from the progenitor star, on the one hand, to the relic compact object, expelled gas, and  ~0.2 M of mass-equivalent neutrino loss, on the other. The challenge, however, is connected to explaining the exact physical reason of the necessary asymmetry.

Although only one percent anisotropy of the several 1053 erg of gravitational binding energy released through neutrinos in the NS formation event would be sufficient to account for an NS recoil velocity of about 300 km s-1, it turns out that an asymmetry of this size or greater is extremely difficult to obtain. Superstrong (≳1016 G) internal, ordered magnetic fields, i.e., fields with a very strong dipolar component, would be needed in the NS. Such fields could impose a direction dependence on the matter-interactions of active neutrinos in the dense NS core (e.g., Bisnovatyi-Kogan 1996; Lai & Qian 1998; Sagert & Schaffner-Bielich 2008), or they could create an asymmetry of resonant Mikheev-Smirnov-Wolfenstein flavor transformations between active neutrinos and sterile neutrinos with masses of several keV (e.g., Fryer & Kusenko 2006; Kusenko 2009; Kishimoto 2011). Also a radiatively driven magnetoacoustic instability, termed neutrino-bubble instability, has been proposed as a possible cause of sizable neutrino emission anisotropies (Socrates et al. 2005). However, all these scenarios require the favorable combination of several unestablished ingredients.

In contrast, NS kick theories based on anisotropic mass ejection during the explosion can be motivated by the observational fact that core-collapse SNe and their gaseous remnants exhibit large-scale asphericity essentially in all observed cases (Leonard et al. 2006; Wheeler & Akiyama 2010; Vink 2012). While the fundamental requirement of the “hydrodynamical” recoil mechanism is thus well consolidated, the exact origin and size of the mass ejection asymmetry needs to be worked out and requires hydrodynamical modeling for quantitative answers. Janka & Müller (1994), using numerical simulations and analytic estimates, showed that the interaction of the proto-neutron star (PNS) with the violent mass flows in the neutrino-heated postshock layer prior to and around the onset of the SN explosion cannot explain NS velocities of many 100 km s-1. This conclusion is supported by more recent three-dimensional (3D) simulations by Fryer & Young (2007). These works relied on asymmetries that developed only after core bounce by neutrino-driven convection and the standing-accretion shock instability (SASI; Blondin et al. 2003; Foglizzo 2002). These hydrodynamic instabilities can grow from random seed perturbations, which can be expected to be present with small amplitudes in any convectively stirred stellar environment. In contrast, Burrows & Hayes (1996) assumed the existence of a large, dipolar asymmetry already in the precollapse stellar core as the relic of convective shell-burning during the late evolution stages of the progenitor star. With the density accordingly reduced in a wide polar funnel, they obtained a pronounced pole-to-pole asymmetry of the explosion, because the ejecta could accelerate more readily to higher velocities on the side of the low-density funnel. As a consequence, the NS received a kick of more than 500 km s-1 opposite to the direction of the stronger explosion. Although the assumed density reduction might be motivated by the observation of vigorous dynamical interaction between different burning shells in combination with large (low-order mode) deviations from spherical symmetry in two-dimensional (2D) stellar calculations for a stage close to core collapse (Arnett & Meakin 2011), the amplitude and structure of such density and velocity perturbations in precollapse stellar cores are still highly uncertain. The determination of these perturbations requires full-sphere 3D simulations over much longer periods of precollapse evolution than the existing 2D models.

Scheck et al. (2004, 2006) therefore also started with minimal and least predetermining assumptions about inhomogeneities in the progenitor before collapse. Like Janka & Müller (1994) they applied small random perturbations (typically 0.1−1% of the density and/or velocity) to trigger the development of large-scale, nonradial flows in the postshock accretion layer, but in contrast to the previous work they could follow the asymmetric, neutrino-powered SN blast wave well beyond the onset of the explosion for evolution times of at least one second (corresponding to shock radii of  ~10 000 km and more). Their results for a large set of 2D simulations demonstrated that due to the action of neutrino-driven convection and low-mode SASI activity, the postshock shell can be expelled with sufficiently large asymmetry to accelerate the NS over a timescale of seconds to velocities up to  ~1000 km s-1 and more. The NS recoil develops in the direction opposite to the side of the strongest explosion. It is mediated mainly by the gravitational forces between the compact remnant and the most slowly moving, massive ejecta “clumps”, whose radial propagation lags behind the rest of the expelled postshock material. The long-distance coupling by gravitational forces can still act on the NS long after hydrodynamical interactions have ceased. The influence of long-range gravity is therefore much more efficient in accelerating the NS into a well defined direction than the random kicks exerted by convective downdrafts during the accretion phase. The impacts of accretion flows buffet the accretor in varying directions. Such a random-walk like process is ineffective in producing high recoil velocities. Recent 2D results of Nordhaus et al. (2010, 2012) lend support to the “gravitational tug-boat mechanism”1, and Wongwathanarat et al. (2010b) presented a first, small set of 3D simulations which showed that NS recoil velocities of at least  ~500 km s-1 can also be obtained in the absence of the artificial contraint to axisymmetry associated with 2D modeling.

In this paper we discuss results of a significantly extended set of 3D simulations for nonrotating 15 M and 20 M progenitor stars, some of which were continued for nearly 3.5 s beyond bounce to follow the saturation of the NS kick velocity. As in the previous works by Scheck et al. (2004, 2006) and Wongwathanarat et al. (2010b), we initiate the neutrino-driven explosions artificially by imposing suitably chosen neutrino luminosities at the inner grid boundary. The latter is placed well inside the neutrinosphere (at neutrino optically depths between  ~10 and some 100), and neutrino effects on the computational grid are treated by an approximate description of the neutrino transport. We also present estimates of the spins that are acquired by the NS. Moreover, we argue and demonstrate by results that the explosive production of heavy elements during the first seconds of the SN blast becomes highly anisotropic in the case of large NS kicks. In particular nuclei between silicon and the iron group are ejected in considerably larger amounts in the hemisphere pointing away from the NS velocity. We therefore predict that the bulk of such nucleosynthetic products should possess a momentum opposite to that of a high-velocity NS. Confirming such a hemispheric momentum asymmetry in SN remnants with determined NS velocities would provide extremely valuable information about the underlying NS kick mechanism. In particular, it might allow one to discriminate our gravitational tug-boat mechanism from scenarios like the neutrino-driven kicks discussed by Fryer & Kusenko (2006), who expect the stronger explosion (and thus enhanced explosive nucleosynthesis) in the direction of the NS motion. High-resolution spectral investigations of young SN remnants like Puppis A (Katsuda et al. 2008, 2010), G11.2-0.3 (Moon et al. 2009), and Cassiopeia A (Isensee et al. 2010; DeLaney et al. 2010; Rest et al. 2011; Hwang & Laming 2012) in different wavebands may offer a promising perspective.

The outline of our paper is as follows. In Sect. 2 we describe our numerical setup, employed methods, and the computed models. In Sect. 3 we present our results for NS kicks, discuss their connection to blast-wave asymmetries (Sect. 3.2), give a detailed description of the physics of the NS acceleration mechanism (Sect. 3.3), provide analytic estimates of the achievable NS velocities (Sect. 3.4), and address possible progenitor dependences (Sect. 3.5) and the contribution from anisotropic neutrino emission (Sect. 3.6). In Sect. 4 we evaluate our models concerning the spin-up of the nascent NSs, and in Sect. 5 we discuss our expectations for the connection between NS kicks and observable asymmetries of the explosive nucleosynthesis of heavy elements. Finally we summarize and draw conclusions in Sect. 6.

thumbnail Fig. 1

Axis-free Yin-Yang grid plotted on a spherical surface demonstrating symmetric geometry of the Yin (red) and Yang (blue) regions. Both grids are the low-latitude part of a spherical polar grid, and therefore contain grid cells which are almost equidistant in angular directions.

Open with DEXTER

2. Numerical methods and setup

2.1. Code, grid, and boundary conditions

Our numerical code is based on the explicit finite-volume, Eulerian, multi-fluid hydrodynamics code Prometheus (Fryxell et al. 1991; Müller et al. 1991a,b). It integrates the multi-dimensional hydrodynamic equations using dimensional splitting following Strang (1968), piecewise parabolic reconstruction (PPM; Colella & Woodward 1984), and a Riemann solver for real gases (Colella & Glaz 1985). Inside grid cells with strong grid-aligned shocks, fluxes computed from the Riemann solver are replaced by the AUSM+ fluxes of Liou (1996) in order to prevent odd-even decoupling (Quirk 1994). The code treats advection of nuclear species by employing the consistent multi-fluid advection (CMA) scheme of Plewa & Müller (1999).

The code employs an axis-free overlapping “Yin-Yang” grid (Kageyama & Sato 2004) in spherical polar coordinates, which was recently implemented into Prometheus, for spatial discretization (Wongwathanarat et al. 2010a). Figure 1 shows the Yin-Yang grid plotted on a spherical surface. The Yin-Yang grid helps relaxing the restrictive CFL-timestep condition and avoiding numerical artifacts near the polar axis. Each grid patch is simply a part of the usual spherical polar grid, and is geometrically identical. This allows us to still make use of the “ray-by-ray” neutrino transport scheme already available in our code. On the other hand, it also means that the Yin-Yang grid still retains the problem of the singularity at the coordinate origin that is present in a spherical polar grid. This, however, is irrelevant in our simulations from the technical point of view because, as in Scheck et al. (2006), we excise the central region of the computational domain interior to a chosen value of the time-dependent inner grid boundary Rib. Our standard grid configuration consists of 400(r)  ×  47(θ)  ×  137(φ)  ×  2 grid cells, corresponding to an angular resolution of 2° and covering the full 4π solid angle. The grid configuration with reduced angular resolution (distinguishable from the 2° configuration by the suffix “lr” in the model name) consists of 400(r)  ×  20(θ)  ×  56(φ)  ×  2 zones, corresponding to an angular resolution of 5°. The radial grid has a constant spacing of 0.3 km from the inner grid boundary Rib up to a radius Rc of 80 km for the progenitors W15, N20, and B15, and 115 km for the progenitor L15 (see Sect. 2.5). Beyond this radius the radial grid is logarithmically spaced. The outer grid boundary Rob is at 18 000–22 000 km, which is sufficient to prevent the supernova (SN) shock from leaving the computational domain during the simulated time. Hydrostatic equilibrium is assumed at the inner grid boundary Rib, while a free outflow boundary condition is employed at the outer one.

Table 1

Chosen values of parameters for all simulated models (see Sect. 2.5 for definitions).

2.2. Input physics

Self-gravity is taken into account by solving Poisson’s equation in its integral form, using an expansion into spherical harmonics as in Müller & Steinmetz (1995). The monopole term of the potential is corrected for general relativistic effects as described in Scheck et al. (2006) and Arcones et al. (2007). The inner core of the PNS with densities well above those of the neutrinospheric layer is excised and replaced by a point mass pinned to the coordinate origin. Thus the PNS core is not allowed to move. The implication of this fact will be discussed in Sect. 3.3. The cooling of the PNS is then described by neutrino emission properties (luminosities and mean spectral energies) that are prescribed as boundary conditions at Rib. The boundary neutrino luminosities are assumed to be constant during the first second of the evolution and decline with time as t−3/2 afterwards (cf. Scheck et al. 2006; Arcones et al. 2007). The mean energies are prescribed as functions of the temperature in the first radial cell of the computational grid as in Ugliano et al. (2012). The contraction of the PNS is mimicked by a radial grid movement with defined velocity according to Scheck et al. (2006); detailed information on the chosen parameters can be found in Table 1. There is no mass flow through the inner grid boundary as long as we follow the evolution of the PNS surface layers.

“Ray-by-ray” neutrino transport and neutrino-matter interactions are approximated as in Scheck et al. (2006) by radial integration of the one-dimensional (spherical), grey transport equation for νe, , and heavy-lepton neutrinos and all angular grid directions (θ, φ) independently. This approach allows us to take into account angular variations in the neutrino fluxes produced by the matter accreted onto the newly forming NS. The accretion luminosity adds to the neutrino fluxes imposed spherically symmetrically at the inner grid boundary so that the outgoing luminosities are considerably different from the values assumed at the inner grid boundary (see, e.g., Scheck et al. 2006; Arcones et al. 2007). The neutrino spectra are assumed to have Fermi-Dirac shape with chemical potentials that are equal to the equilibrium values in optically thick regions and constant, limiting values in the free streaming regime. The spectral temperatures are computed from inverting the ratio of neutrino energy and number fluxes, whose transport is solved simultaneously. Details can be found in Scheck et al. (2006). The tabulated equation of state (EoS) of Janka & Müller (1996) is used to describe the stellar fluid. It includes arbitrarily degenerate and arbitrarily relativistic electrons and positrons, photons, and four predefined nuclear species (n, p, α, and a representative Fe-group nucleus) in nuclear statistical equilibrium.

2.3. Nucleosynthesis

In order to follow approximately the explosive nucleosynthesis, we solve a small α-chain reaction network, similar to the network described in Kifonidis et al. (2003). For our simulations we choose a subset of nine α nuclei, 4He, 12C, 16O, 20Ne, 24Mg, 28Si, 40Ca, 44Ti, 56Ni, and an additional “tracer nucleus”. We omit 32S, 36Ar, 48Cr, and 52Fe from the complete α-series in order to reduce the network size and thus to save computing time. The two gaps between 28Si–40Ca and 44Ti–56Ni are bridged by combining the intermediate reaction steps and estimating the new effective reaction rates by taking the slowest rate of the intermediate reactions. The tracer nucleus is produced via the reaction 44Ti(3α,γ)56Ni within grid cells whose electron fraction Ye is below 0.49. Such conditions are found in the neutrino-heated ejecta of our models. The tracer nucleus represents iron-group species that are formed under conditions of neutron excess. It thus allows us to keep track of nucleosynthesis in neutron-rich regions. However, it should be noted that some fraction of the tracer material may actually be 56Ni, because our approximations in the neutrino transport tend to underestimate Ye in the neutrino-heated ejecta, while more sophisticated energy-dependent neutrino transport yields (slightly) proton-rich conditions in the neutrino-processed material expelled during the early explosion (e.g., Pruet et al. 2006; Fröhlich et al. 2006; Fischer et al. 2010; Müller et al. 2012). The network is solved in grid cells whose temperature is within the range of 108 K and 8 × 109 K. We neglect the feedback from the network composition to the EoS and the effect of the nuclear energy release on the hydrodynamic flow. The energy release by nuclear reactions is of minor relevance for the dynamics because the production of 0.1 M of 56Ni means a contribution of only  ~1050 erg to the explosion energy. It is important to note that our model parameters are chosen to give energetic explosions already by neutrino energy input. Not using an exact nuclear composition in the EoS is an acceptable approximation because at neutrino-heated ejecta conditions the contributions of nuclei to pressure, energy density, and entropy are dwarfed by those of electrons, positrons, and photons. To prevent too small nuclear burning timesteps we do not perform any network calculations above 8 × 109 K. When freeze-out from nuclear statistical equilibrium (NSE) is followed in ejecta cooling down from above 8 × 109 K, we start from pure α-composition, which is compatible with the high-temperature NSE composition in the absence of free neutrons and protons, in particular since Ye = ne/nB has values close to 0.5 in all of the ejecta. As the hot ejecta expand and their temperature decreases these α-particles recombine to produce heavier α-nuclei considered in our α-network, among them 56Ni. Once the temperature drops below 108 K, all nuclear reactions are switched off because they become too slow to change the nuclear composition on the explosion timescale. Nuclear burning in shock-heated matter is also described by our α-network.

2.4. Long-time simulations

To follow the propagation of the SN shock wave until very late times, we map our results after 1.1–1.4 s onto a new computational grid, whose inner and outer radial boundaries are placed at 500 km and 3.3 × 108 km, respectively. The latter value is near the stellar surface. At the inner grid boundary we assume the inflow of a spherically symmetric neutrino-driven wind, which is considered to be the consequence of ongoing neutrino emission and corresponding energy deposition near the PNS. We do not directly simulate neutrino-matter interactions during our long-time runs, because details of the neutrino physics become less important during the late evolution phases when the explosion conditions are already determined. By including a wind inflow, however, we continue to account for the consequences of neutrino heating at the PNS surface. Since the imposed wind has spherical symmetry it will not affect the anisotropy of the leading ejecta that we intend to investigate. The hydrodynamic inflow quantities are obtained as the angular averages of the flow properties at 500 km at the time of the mapping and kept constant for another 2 s of evolution.

The wind is replaced by a free outflow boundary condition afterwards. At the outer boundary we apply a free outflow condition at all times. Since we keep the angular grid resolution unchanged but increase the number of radial zones, our computational grid grows to 1200(r)  ×  47(θ)  ×  137(φ)  ×  2 zones. The radial resolution is approximately Δr/r ≈ 0.01. We successively discard the innermost radial zone and thus move the inner grid boundary outward to a new location whenever its radius becomes less than 2% of the minimum shock radius. Moving the inner grid boundary to a larger radius helps relaxing the CFL timestep and allows us to follow the evolution to very late times with acceptable computational costs.

Also during the long-time runs we consider nucleosynthesis processes and account for the self-gravity of the stellar gas by a Newtonian description. We switch to the tabulated EoS of Timmes & Swesty (2000), appropriate for the lower density and temperature values present in the outer stellar layers (not covered by the EoS table of Janka & Müller 1996). Moreover, besides considering the contributions from arbitrarily degenerate and relativistic electrons and positrons and from a photon gas, it treats the nucleonic gas components by a mix of ideal Boltzmann gases with 11 nuclear species: p, 4He, 12C, 16O, 20Ne, 24Mg, 28Si, 40Ca, 44Ti, 56Ni, and X. This is sufficient for a reasonably good representation of the composition in the outer core, mantle, and envelope of the progenitor.

Table 2

Explosion and NS properties for all models (see Sect. 3.1 for definitions).

2.5. Computed models

We investigate three 15 M models and one 20 M progenitor denoted as W15, L15, B15, and N20, respectively. W15 is based on the nonrotating 15 M model s15s7b2 of Woosley & Weaver (1995), L15 corresponds to a star evolved by Limongi et al. (2000), B15 is a blue supergiant (SN 1987A) progenitor of Woosley et al. (1988), and N20 was computed up to core collapse by Shigeyama & Nomoto (1990). The 15 M progenitors W and L were followed through collapse to 15 ms after bounce with the Prometheus-Vertex code (Rampp & Janka 2002) in one dimension (Marek & Buras, priv. comm.). The Prometheus-Vertex code is an Eulerian multi-D hydro-code (based on a finite volume method) coupled to a two-moment closure scheme with variable Eddington factor for neutrino transport. Using the same code, the progenitor N20 was evolved until 11 ms after bounce (Marek, priv. comm.), whereas the model B15 was collapsed in a 1D simulation (Bruenn 1993) using a Lagrangian hydro-code with multi-group flux-limited neutrino diffusion (Bruenn 1985). Both simulation codes are Newtonian, but Prometheus-Vertex employs correction of the gravitational potential due to general relativistic effects (Marek et al. 2006). To break spherical symmetry, random seed perturbations with an amplitude of 0.1% (and cell-to-cell variations) are imposed on the radial velocity (vr) field2. Explosions with chosen energy are initiated by neutrino heating at a rate that depends on suitable values of the neutrino luminosities assumed at the lower boundary.

We computed in total 20 models varying the neutrino luminosities imposed at the inner grid boundary and thus the resulting explosion energies (Table 2). The choice of parameters for each model is listed in Table 1. These are Lνe,ib and , the luminosities of electron neutrinos and electron antineutrinos at the inner grid boundary, respectively3, Rib, the initial radius of the inner grid boundary, , the asymptotic radius of the inner grid boundary as t → ∞, tib, the timescale for the boundary contraction, Rob, the radius of the outer grid boundary, and Rc, the radius up to where the radial grid resolution is kept constant at 0.3 km. Some of our models (W15-1 and W15-2; B15-1 and B15-2) differ only by the initial seed perturbations. The evolution including our treatment of neutrino transport is followed until 1.3 s after bounce for the W15 and N20 models, while the L15 simulations are continued until 1.4 s postbounce and the B15 runs until about 1.1 s after bounce. The time at which we stop our calculations differs from progenitor to progenitor depending on how fast the explosion sets in. Explosions in the L15 models are the most delayed ones, while the B15 models show the fastest explosion. The simulations with detailed neutrino treatment are stopped when the maximum SN shock’s radius is close to Rob. Models which are simulated with the reduced angular grid resolution of 5° are indicated by the suffix “lr” appended to their model names. In addition, a subset of 9 models is continued into the long-time evolution as described in Sect. 2.4 in order to study the further acceleration of the PNS and the explosive nucleosynthesis behind the outgoing shock.

thumbnail Fig. 2

Entropy-isosurfaces (left) of the SN shock (grey) and of the high-entropy bubbles (green), and entropy distribution in a cross-sectional plane (right) at t = 140,248,515 ms, and 1.3 s for model W15-2. The viewing direction is normal to the xz-plane. Length scales are given by the yardsticks below each image. Small Rayleigh-Taylor mushrooms start growing at  ~100 ms (top panels). The rising high-entropy bubbles merge and rearrange to form larger bubbles by the time the explosion sets in (second from top). The NS starts accelerating due to the asymmetric distribution of the ejecta, the acceleration reaching its maximum at  ~500 ms (third from top). At this epoch, the ejecta show a clear dipolar distribution with more dense, low-entropy material concentrated in the hemisphere containing the kick direction. The final NS kick direction has already found its orientation, and its projection onto the xz-plane is shown by the white arrows pointing to the lower right direction. The bottom panels show the entropy structure when the essentially spherically symmetric neutrino-driven wind has developed, visible as the green central region enclosed by the highly aspherical wind-termination shock (lower right panel).

Open with DEXTER

thumbnail Fig. 3

Entropy-isosurfaces (left) of the SN shock (grey) and of the high-entropy bubbles (green), ray-casting images of the density (middle), and entropy distribution in a cross-sectional plane (right) for models W15-1 (top) and L15-2 (bottom) at about 1.3 s and 1.4 s after core bounce. The viewing directions are normal to the plane spanned by the NS kick and spin vectors of each model, which also defines the plane for the entropy slices. The SN shock has an average radius of 13 000 km and 14 000 km for models W15-1 and L15-2, respectively (a length of 5000 km is given by yardsticks below the left images). The kick and spin directions are indicated by the white and black arrows, respectively, in the middle figures. The NS locations are marked by black crosses in the right plots. The images in the middle correspond roughly to the projections of the density distribution on the observational planes. Dilute bubble regions are light-colored in white and yellow, while dense clumps appear more intense in reddish and bluish hues. The purple circular areas around the NS represent the dense inner region of the essentially spherically symmetric neutrino-driven wind. The wind is visible in green in the right images and is bounded by the aspherical wind-termination shock. The NS is clearly displaced from the geometrical center of the expanding shock in the direction of the kick vector pointing to the lower left for model W15-1. The NS in model L15-2 has a much smaller kick velocity than in model W15-1, and does not show any clear displacement but remains located roughly at the center of the expanding ejecta shell.

Open with DEXTER

3. Neutron star kicks

Wongwathanarat et al. (2010b) presented results of NS kicks for a small set of 3D computations conducted until 1.3–1.4 s after bounce. Here, we describe the evolution of these models for a much longer time as well as more results of a significantly enlarged model set that includes runs for three 15 M progenitors and the 20 M star. Moreover, we analyze our numerical results in the light of theoretical considerations and simple analytical estimates (extending the brief discussion in Wongwathanarat et al. 2010b, and sharpening arguments made by Scheck et al. 2006; and Nordhaus et al. 2010). We also evaluate our simulations for recoil velocities caused by anisotropic neutrino emission in comparison to the far dominant kicks by asymmetric mass ejection.

3.1. Simulation results

Figure 2 presents snapshots of model W15-2 at four different epochs: 0.140, 0.248, 0.515, and 1.30 s after bounce. Violent buoyancy and mass overturn is triggered by neutrino heating around the NS. At roughly 100 ms, Rayleigh-Taylor mushrooms begin to grow from the imposed seed perturbations in the convectively unstable layer that builds up between the gain radius (that separates neutrino-cooling and -heating regions) and the stalled SN shock. These high-entropy bubbles evolve in a highly non-stationary manner, start rising, growing, merging, partially collapse again and emerge once more to inflate to larger sizes. Along with the slowly expanding shock their angular diameter increases, adjusting roughly to the radial extension of the convective layer. Aided by convective overturn and global shock motions, possibly associated with activity due to the SASI, the delayed, neutrino-driven explosion sets in at about 250 ms postbounce.

To identify the explosion epoch, we perform a time-dependent evaluation of the explosion energy Eexp defined as the sum of the total (i.e., internal plus kinetic plus gravitational) energy of all grid cells where this energy is positive. We define the explosion time texp as the moment when the computed explosion energy exceeds 1048 erg. This roughly coincides with the time when the average shock radius exceeds 400–500 km.

This is also the time when the NS develops a growing kick velocity along a well defined direction. While the kick direction was randomly fluctuating before, the acceleration vector now settles and the NS velocity grows (see Wongwathanarat et al. 2010b). A cross-sectional view of the entropy distribution clearly shows a dipolar asymmetry with dense, low-entropy material concentrated at small radii primarily in the hemisphere to which the kick vector is pointing (Fig. 2, third row and also Fig. 3, upper row). At a later stage, the pattern of the nonradial structures freezes and the ejecta continue to expand in a nearly self-similar manner. An entropy slice shows the emergence of an essentially spherically symmetric neutrino-driven wind enclosed by an aspherical wind-termination shock. On the side where the explosion is stronger, which manifests itself by faster shock expansion and higher entropy behind the SN shock (see the 10−11 o’clock position in the bottom right panel of Fig. 2 and the upper right panel of Fig. 3), the wind is shocked to higher entropies because it passes the termination shock at larger radii and thus with higher velocity (see also Arcones et al. 2007 and Arcones & Janka 2011).

Table 2 summarizes explosion and NS properties for all computed models: Mns, texp, Eexp, vns, ans, vns, αkν, Jns, αsk, and Tspin denote the NS mass, time of the onset of the SN explosion, explosion energy, NS kick velocity by the gravitational tug-boat mechanism, corresponding NS acceleration, additional neutrino-induced kick velocity, angle between gravitational and neutrino-induced NS kick, NS angular momentum, angle between NS spin and kick directions, and estimated final NS spin period. While results for all models are provided after 1.1–1.4 s of post-bounce evolution, and give the NS velocity and acceleration at 3.1–3.4 s after bounce for the subset of the simulations that was continued until this later time.

The NS mass Mns is defined as the baryonic mass enclosed by the NS radius Rns, which is considered as the angular dependent radial location where the density is 1011 g cm-3. In estimating the NS recoil velocity, we assume conservation of the linear momentum of the whole progenitor star. Since the high-density core of the PNS is excised in our simulations, the NS cannot move when it exchanges momentum with the surrounding matter. Thus it behaves like an object with infinite inertial mass, absorbing momentum in analogy to a wall that reflects a bouncing ball. Nevertheless, following Scheck et al. (2006), we can estimate the NS velocity vns from the negative of the total linear momentum of the ejecta gas as (1)where is the integral of the gas momentum over the computational volume between PNS radius and outer grid boundary4.

In our simulations we find NS kick velocities ranging from about 30 km s-1 (L15-3) to roughly 440 km s-1 (W15-6) at 1.1−1.4 s after bounce. The NS acceleration ans is obtained by numerically differentiating the NS velocity calculated from Eq. (1) and listed in Table 2. Because of the nonvanishing acceleration, the recoil velocity at 1.1–1.4 s p.b. is usually not the final NS velocity. Only models ending with a low kick (like L15-2 and L15-3) may have reached nearly the asymptotic value already after this period of time. In the three long-time runs W15-1, W15-2, and W15-6, we find a significant increase in the recoil velocity between 1.3 and 3.3 s, and even then the acceleration continues, though at a considerably smaller rate (compare the values of vns and ans with and without the superscript “long” in Table 2). Although we evolved several of the long-time models to much later times we do not compute the NS velocities beyond 3.4 s, because the subsequent shift of the inner grid boundary and the application of a free outflow condition, can lead to anisotropic momentum flow across the boundary can occur. The latter would have to be accounted for with high precision if the argument of momentum conservation were to be used for estimating the NS recoil also later on.

thumbnail Fig. 4

Long-time evolution of the NS velocity (top) and of the corresponding NS displacement (bottom) for models W15-6, W15-1, W15-2, L15-1, and L15-2. The displacement of the kicked NS from the coordinate center, although not taken into account in our simulations, can be estimated by time integration of the NS velocity. Within 3.3 s the NS in model W15-6 exceeds a velocity of 700 km s-1 and continues to be accelerated at a modest rate. Note also that the NS velocities in models L15-1 and L15-2 have nearly reached their terminal values after 3 s. In all cases the NS shift relative to the coordinate center is negligibly small before the explosion takes off, and even in model W15-6 with its most extreme NS kick it is always dwarfed by the diameter of the shock wave and postshock ejecta, which expand with typically  ~10 000 km s-1.

Open with DEXTER

In Fig. 3 the asymmetry of the SN shock, the structure of the high-entropy bubbles, and the inhomogeneous, dense ejecta shell behind the expanding shock are compared for two models at 1.3–1.4 s. One model, W15-1 (top), has a large kick, the other model, L15-2 (bottom), a small one. The analogous combination of plots was provided for model W15-2, our second highest kick case, in Fig. 3 of Wongwathanarat et al. (2010b). The three panels display different aspects of the anisotropic mass distribution. While in the left image the geometry is visible in overview by three-dimensional entropy-isosurfaces, the middle panel gives an impression of the asymmetry of the gas (column) density visualized by a ray-casting technique. The viewing direction is perpendicular to the plane defined by the NS kick (white) and spin (black) vectors. The right plot shows the entropy distribution in that plane. The ray-casting images visualize the density projected along the lines of sight. They show that the kick vector points to the region with the highest concentration of dense clumps. Moreover, the NS in model W15-1 is visibly displaced from the geometrical center of the ejecta and of the SN shock in the recoil direction; this effect is even more pronounced for W15-2. In contrast, in the low-kick case of L15-2 the high-entropy bubbles are of relatively similar sizes on all sides. The expansion of the SN shock is therefore nearly equally strong in all directions, i.e. the shock has no strong dipolar deformation (though a quadrupolar asymmetry mode is visible), and the NS is close to the center of the ejecta shell. Correspondingly, the wind-termination shock is relatively spherical in this model. While there are denser, clumpy structures near the one o’clock position of model L15-2, there are also similarly massive, though less concentrated filaments in the opposite hemisphere. The explosion asymmetry has a dipolar component that is not large enough to produce a great pulsar kick.

thumbnail Fig. 5

Time evolution of the normalized spherical harmonic expansion coefficients for different (ℓ,m) modes computed for the shock surface Rs(θ,φ) (left) and for the surface mass density D(θ,φ) (Eq. (2); right) of models W15-2 (top) and L15-2 (bottom). The polar axis of the coordinate system chosen for the expansion is aligned with the final NS kick direction. Note the dominance of the (ℓ,m) = (1,0) mode for model W15-2 with the highest kick velocity, while this dominance is absent in the case of model L15-2.

Open with DEXTER

Figure 4 displays the growth of the NS velocities in five of our nine long-time runs over the whole period of 3.3–3.4 s. In the most extreme cases, models W15-6 and W15-2, the NS velocity reaches an asymptotic value over 700 km s-1 and near 600 km s-1, respectively. Several of the W15-models therefore develop kicks clearly higher than the peak of the measured velocity distribution of young pulsars. Only more simulations for a wider variety of progenitors, however, will be able to provide answers to the questions whether velocities of even 1000 km s-1 and more can be obtained by the discussed gravitational tug-boat mechanism, and whether and how such extreme recoils might be connected to particular features of the progenitor stars or especially favorable conditions of the explosion itself. Analytic estimates in Sect. 3.4 will show that the mechanism is likely to be viable of producing kicks well in excess of 1000 km s-1 for reasonable assumptions about the explosion asymmetry. Yet it is not clear that the magnitude of the kick is merely a stochastic result of the chaotic growth of asymmetries by hydrodynamic instabilities in the SN core. Some trends of our results seem to indicate that certain properties of the progenitor and the explosion could be supportive for large kicks (see also Sect. 3.5).

3.2. Mode analysis

The causal connection between the explosion asymmetries, especially the density asymmetries in the ejecta shell, and the natal NS kicks, which is visually suggested by Figs. 2 and 3 as well as by Fig. 3 of Wongwathanarat et al. (2010b), can be quantitatively supported by an expansion of variables that characterize the mass distribution in the postshock layer into spherical harmonics.

Testing different possibilities, we ascertained that the most suitable quantity for this purpose is the “surface mass density” (2)where Rs(θ,φ) is the direction-dependent radius of the SN shock. The integrand of Eq. (2) means that mass elements Δm at different distances from the NS are weighted according to Δm/r2, reflecting the importance of their contribution to the NS acceleration by gravitational effects. In contrast, the shock radius Rs(θ,φ) or the integrated mass per solid angle, (3)put too much weight on the region just behind the shock, which is not necessarily the most relevant one for accelerating the NS. The NS kick is mainly determined by the density inhomogeneities in the deeper ejecta regions. Although a large kick velocity is typically also linked to a pronounced  = 1 displacement of the shock surface relative to the grid center (in our most extreme cases the dipole amplitude can reach nearly 10% of the monopole one, see Fig. 5, left panel), we find that the inverse conclusion does not always hold: A dominant  = 1 shock asphericity does not necessarily imply a strong NS recoil. In fact, we observe the  = 1 spherical harmonics of the shock surface to be larger than all higher-order multipole amplitudes at late postbounce times for the majority of our models (except for the few cases like L15-2, L15-3 or N20-2, all of which are among our runs with the lowest NS kicks). For a discriminating diagnostic quantity we therefore prefer to evaluate D(θ,φ), whose mode amplitudes are more conclusive about ejecta asymmetries that are causal for large recoil velocities of the compact remnant. Large kicks correlate with an  = 1 mode of D(θ,φ) that clearly dominates the higher modes during (at least) a longer period of the postbounce evolution, whereas for cases with small kicks the  = 1 mode is not much stronger at any time than the other modes.

The decomposition of D(θ,φ) in spherical harmonics is written as (4)where are the expansion coefficients and the spherical harmonics are functions of the associated Legendre polynomials  as (5)with (6)Multiplying Eq. (4) by the complex conjugate of the spherical harmonic, , and integrating over the solid angle, the expansion coefficient is found to be (7)We choose the coordinate system for the expansion such that the polar axis is aligned with the final NS kick direction. Before performing the integration of Eq. (7), we interpolate, using bi-linear interpolation, the surface mass density computed on the Yin-Yang grid onto the new coordinate system. Figure 5, right panel displays the resulting expansion coefficients of different (ℓ,m) modes versus time for models W15-2 and L15-2. In the former, high-kick model the amplitude of the (ℓ,m) = (1,0) mode reaches its maximum at about 400 ms after bounce, i.e., approximately at the same time when the NS acceleration becomes maximal (Fig. 6). After this time, at t ≳ 0.4 s postbounce, the (ℓ,m) = (1,0) mode clearly dominates. Such a behavior is absent in the low-kick model L15-2.

Consistent with what is visible in Figs. 2, 3 and Fig. 3 of Wongwathanarat et al. (2010b), the mode analysis shows that a dominant  = 1 (dipolar) asymmetry of the mass distribution in the ejecta correlates with a high pulsar kick along the dipole axis. In model L15-2 all of the  = 1 and  = 2 modes displayed in the lower right panel of Fig. 5 exhibit similar amplitudes and thus the gravitational pull of the ejecta does not develop a preferred direction. The acceleration of the NS remains correspondingly low. This result provides strong support to the described gravitational tug-boat mechanism for NS kicks, which is linked to the anisotropic ejection of matter in the first second of the supernova explosion and which is mediated to the NS on significantly longer timescales by the asymmetric gravitational pull of the expanding ejecta gas. In the next section we will describe in more detail how the asymmetries are established and how the acceleration mechanism works. Moreover, we will present simple toy models to approximately estimate its possible strength in dependence of the properties of the supernova ejecta.

thumbnail Fig. 6

Left: time evolution of the NS acceleration in model W15-2. The acceleration computed by the assumption of momentum conservation, i.e., by the negative derivative of the total linear momentum of the gas surrounding the NS (red; denoted by “mom.consv.”), is compared with the value obtained by the integration of all hydrodynamical and gravitational forces that act on a sphere of radius 1.3Rns around the PNS (black; denoted as “total”). Also shown are the absolute values of the individual contributions of all terms in Eq. (8), namely the momentum transfers by downflows (light blue), outflows (green), anisotropic pressure forces (purple), and the gravitational attraction by the asymmetric ejecta (blue). The last contribution clearly dominate the evolution over long periods. Right: NS velocity components in all coordinate directions and total value of the NS velocity. The two lines in each panel, which correspond to our two methods of computing the NS recoil acceleration in the left panel, can hardly be distinguished from each other.

Open with DEXTER

3.3. Kick mechanism: theory and tests

The NS acceleration that is associated with the asymmetric ejecta distribution developing during the onset of the SN explosion is achieved mainly by gravitational forces in combination with hydrodynamical forces. This will be detailed below. The acceleration proceeds in three steps:

  • (1)

    When violent convective mass flows and possibly SASIsloshing motions conspire to stir the postshock layer, ananisotropy of the mass-energy distribution around the PNSis created. Convective downdrafts, channelling accretionflows from the stalled shock through the neutrino-heating re-gion to the vicinity of the gain radius, get deflected to feed anasymmetric pattern of high-entropy bubbles. The bubblesform, grow, contract again, partly collapse, and reappear in aquasi-chaotic way to become smaller or larger, absorbing lessor more neutrino energy. This stochastic bubble formation,however, does not cause an appreciable recoil of the NS (seeFig. 6 and Wongwathanaratet al. 2010b).

  • (2)

    When the explosion sets in, the shock and postshock gas begin to expand aspherically, driven by the asymmetric inflation of the bubbles. The ejecta gas therefore gains radial momentum and its c.o.m. moves away from the coordinate origin (Figs. 2 and 3, right panels): The ejecta shell acquires a net linear momentum because of different strengths of the explosion in different directions. The initial energy and mass asymmetry is thus transformed into a momentum asymmetry by the conversion of thermal to kinetic energy through hydrodynamical forces. When the expansion timescale becomes shorter than the timescale of lateral mixing, the asymmetric ejecta structures essentially freeze out.

  • (3)

    Because of linear momentum conservation, the NS must receive the negative of the total momentum of the anisotropically expanding ejecta mass. Hydrodynamic pressure forces alone cannot achieve the NS acceleration (Scheck et al. 2006). As long as accretion downdrafts reach the NS, momentum is transferred by asymmetric gas flows. Stronger accretion on the weaker side of the blast and more mass loss in the neutrino-driven wind on the other side cause a recoil opposite to the main explosion direction. However, the largest kick contribution, which continues even after accretion has ended and after the wind outflow has become spherical (Fig. 2, bottom right panel, and Fig. 3, right panels), results from the gravitational pull of the anisotropic shell of ejecta (Scheck et al. 2006; Wongwathanarat et al. 2010b; and also Nordhaus et al. 2010).

Why can our models produce kick velocities in excess of  ~500 km s-1 by hydrodynamically created explosion asymmetries although such a mechanism appeared to be highly disfavored even for the extreme case of dipolar ejecta deformation in an analysis of 2D explosion models by Janka & Müller (1994)? It should be noted that the kick scenario discussed in the latter paper is fundamentally different from the one described here. Janka & Müller (1994) assumed the kicks to originate from an impulsive event, happening during a short period of a few hundred milliseconds around the onset of the anisotropic explosion. Correspondingly, in their estimates they considered the momentum asymmetry of the ejecta shell during this early phase, expecting the kick to be established before the accretion is over. Indeed, also in the simulations described in the present paper the NS velocity at the end of the accretion phase is much lower than 500 km s-1 in all cases. However, in the scenario discussed here the NS acceleration is a long-duration phenomenon, which continues for several seconds even after accretion has ended, because it is promoted by the long-range gravitational forces of the asymmetric, dense ejecta shell that expands away from the NS behind the outgoing SN shock.

To elaborate on this point and to check our estimates of the NS recoil velocity on the basis of Eq. (1), we follow Scheck et al. (2006) and perform an analysis of the different contributions to the NS kick by summing up all forces acting on a central sphere of radius r0 = 1.3Rns (Nordhaus et al. 2010, adopted a similar analysis for their recent 2D results). According to Eq. (5) of Scheck et al. (2006) the time-derivative of the NS momentum can be deduced from the Euler equations as (8)The corresponding NS acceleration is and the NS velocity . The first term on the rhs of Eq. (8) accounts for the variation in the pressure force around the sphere. The second term is the flux of momentum going through the sphere. This term can be decomposed into downflow and outflow contributions depending on the sign of the radial velocity vr being negative or positive, respectively. The third term is the contribution from the long-range gravitational forces exerted on the NS by the anisotropic mass distribution outside the sphere. We note in passing that a change in the radius of evaluation, r0, in a reasonable range around 1.3Rns (i.e., in the neighborhood of the NS surface) leads to different weights of the three summed terms (as noticed recently by Nordhaus et al. 2010) but still yields the same result for the integrated total effect.

Figure 6 compares ans and vns for model W15-2 as calculated by the two different methods, i.e. from the requirement of total momentum conservation (Eq. (1)) and from the integration of the responsible hydrodynamical and gravitational forces according to Eq. (8). NS acceleration as well as velocity exhibit excellent agreement of the results obtained in the two different ways. Between  ~130 and  ~600 ms after bounce the forces by the anisotropic pressure distribution and by the gravitational pull both contribute to the total acceleration on significant levels, while downflows and outflows transfer sizable amounts of momentum mostly during transient phases of high activity when strong accretion streams impact the nascent NS anisotropically. At later stages, i.e., at t ≳ 400 ms when the asymmetry of the ejecta has been established, and in particular at t ≳ 600 ms when accretion to the PNS has ceased and the neutrino-driven wind has become essentially spherically symmetric (Fig. 2), the gravitational force term starts to dominate and therefore accounts for the long-term increase in the NS kick. A major fraction of the final NS velocity is thus caused by the gravitational attraction of the anisotropic ejecta.

thumbnail Fig. 7

Ejecta morphology of the high-kick models W15-2 (left) and W15-6 (right) at approximately 3.3 s after bounce. The NS (indicated by black spheres, enlarged by a factor of 20 for better visibility) are accelerated by the gravitational tug of the asymmetrically distributed, dense ejecta to velocities of 575 km s-1 and 704 km s-1, respectively (green arrows pointing to the right). The ejecta “clumps” with the highest densities, which are mainly responsible for the gravitational attraction, are visualized by the red-colored surfaces of constant density of 104 g cm-3. The small black dot left of the NS in both images marks the coordinate origin. The NS displacement relative to this location was computed by time integration of the NS velocity. After  ~3.3 s it is only 1260 km in W15-2 and 1440 km in W15-6. The neutrino-heated, asymmetrically distributed inner ejecta are bounded by the semi-transparent, beige surface from outside and the semi-transparent, light blue surface from inside. The former surface is defined by a constant 56Ni mass fraction of 5%, the latter surface conincides with the termination shock of the neutrino-driven wind. The outermost, semi-transparent green surface (cut open towards the observer for better visibility of the inner volume) corresponds to the SN shock wave. The average shock radius is approximately 36 000 km at the displayed time. A yardstick of 10 000 km along with a tripod indicating the coordinate directions are provided in the lower right corner of each plot.

Open with DEXTER

Our analysis, compared to the 2D results of Scheck et al. (2006), reveals that in the 3D case the direct hydrodynamical momentum transfer by gas downflows to and outflows from the NS is much less important than in 2D. While in 2D the toroidal nature of all structures favors the formation of large hemispheric differences with inflow to the NS dominating in one hemisphere and outflow in the other (corresponding to an  = 1 geometry), the 3D downflow and outflow pattern is more often and for longer periods of time characterized by higher multipoles. The momentum transferred to the NS by inflows is thus almost balanced by the momentum carried away from the NS by outflows on smaller angular scales, resulting only in a small contribution to the NS kick velocity. The gravitational tug by the asymmetric ejecta distribution is therefore clearly the dominant mediator of the NS acceleration in a time-integrated sense, although it is not necessarily the only important effect, in particular not during the very early phase of the NS recoil and not in 2D. Sizable directional differences of the explosion strength can occur and can lead to large ejecta inhomogeneities and asymmetries also in 3D (most pronounced in models W15-6 as well as W15-1 and W15-2, see Figs. 2 and 3). Large recoil velocities can therefore be produced despite the fact that the pattern of convective bubbles and downdrafts seems to be characterized by the presence of higher-order multipoles.

Allowing the PNS to move (instead of being attached to the coordinate center because of the excision of the central core in our simulations) is unlikely to change fundamental aspects of our results. This expectation is based on a simple fact visualized in Figs. 4 (bottom panel) and 7. Even for the highest computed kick velocities the displacement of the NS from the coordinate origin during the computed evolution of several seconds (at most  ~1500 km) is very small compared to the radial extension of the explosion shock and asymmetric ejecta distribution (typically several 10 000 km). It is therefore not astonishing that our expectation was recently confirmed by Nordhaus et al. (2010, 2012) with 2D simulations, in which the NS was able to move. The minor influence of the NS motion had been tested before by Scheck et al. (2006), who performed a large set of 2D simulations in which some of the kicks were similarly strong as the 3D cases discussed here. For doing the tests, Scheck et al. (2006) applied a Galilei transformation with velocity −vns to the gas outside of the PNS (see Sects. 2.3, 7.2, Appendix B, and Table A.5 in the cited paper), thus mimicing the effects of an accelerating NS moving through the surrounding medium. Naturally, because of the stochastic and chaotic growth of local and global asymmetries, the simulations developed values of vns that could be different for a fixed progenitor and the same choices of boundary conditions (core neutrino luminosities and PNS contraction) and initial perturbations. Nevertheless, the whole set of models with the NS motion taken into account did not exhibit any systematic differences compared to the standard set without the Galilei transformation applied (cf. Tables A.1 and A.5 in Scheck et al. 2006).

3.4. Analytic estimates of the kick magnitude

What are the maximum kick velocities that can be mediated to the NS by gravitational forces of asymmetrically ejected matter? The maximum velocity of more than 700 km s-1 obtained in our present set of 3D models is no upper limit, and recoil velocities beyond 1000 km s-1 were already found in 2D simulations by Scheck et al. (2006). Simple considerations can yield order-of-magnitude estimates that confirm this possibility. In the following we discuss two different, idealized cases that allow one to integrate the gravitational pull over the whole period of acceleration (see also Wongwathanarat et al. 2010b and alternative, simple estimates by Nordhaus et al. 2010).

In the first case we consider a hemispheric asymmetry of the mass distribution in an expanding, spherical ejecta shell. The shell with an average radius rs is assumed to contain an extra clump-like mass Δm in one hemisphere, while an equal mass deficit exists in the other hemisphere. A mass asymmetry of this kind may be established by asymmetric accretion flows during the violently convective postbounce phase. Downflows may, for example, be stopped and blown out again before they can reach down to the PNS surface, or only some part of the downflows is accreted onto the PNS while the rest is deflected and ejected again after neutrino heating. Let us assume that the shell with its time-dependent radius rs(t) expands with a constant velocity vs, starting at an initial radius ri. If the NS is displaced by a distance s from the center of the shell (s = 0 at t = 0), the gravitational forces lead to a NS velocity vns =  via an acceleration given by (9)Using rs = ri + vst and assuming s ≪ rs at all times, the integration of Eq. (9) from t = 0 to t = ∞ yields: (10)where Δm is normalized by 10-3M, ri by 107 cm, and vs by 5000 km s-1. A 10-3M hemispheric asymmetry of the shell expanding with a constant velocity of 5000 km s-1 can thus pull the NS to a velocity of 540 km s-1 (this result was already mentioned in Wongwathanarat et al. 2010b). Ejecta asymmetries can therefore very effectively mediate a long-lasting acceleration of the NS, which is pulled in the direction of the higher mass concentration. According to Eq. (10) the kick becomes larger for a lower expansion velocity, because the gravitational tug from the asymmetric ejecta shell remains strong for a longer time. While a value of 5000 km s-1 tends to be on the high side for the postshock ejecta shell during the relevant times, a clumpy anisotropy of 10-3M means an asymmetry of the mass distribution of only  ≲1% in a shell that typically contains 10-1M or more. Asymmetries of several percent appear well within reach and therefore velocities of more than 1000 km s-1 as well. It is interesting that the first-order terms of s/rs cancel in Eq. (9). This means that the NS recoil is affected by the NS motion and displacement from the explosion center only on the level of second-order correction terms of s/rs ≪ 1.

In the second case we consider a difference of the expansion velocity of two clumps Δm of the ejecta in both hemispheres instead of a hemispheric mass difference. This means that we assume that these clumps have been accelerated differently strongly and propagate away from the center of the blast according to rs,1 = ri + vs,1t and rs,2 = ri + vs,2t with vs,1 < vs,2. Taking s to be again the displacement of the NS from the blast center, the compact remnant in this situation experiences a gravitational acceleration (11)To lowest order in s/rs ≪ 1, time integration from 0 to ∞leads to (12)where we have introduced the definitions Δvs = vs,2 − vs,1 and . In the last expression we have again normalized Δm by 10-3M, ri by 107 cm, and by 5000 km s-1. The acceleration is opposite to the faster expanding hemisphere, i.e. it is against the direction of the stronger explosion, and of the same order of magnitude as the result of Eq. (10), if we assume a velocity asymmetry of 100%, . This, however, is on the very extreme side, and seems more realistic. According to Eqs. (11) and (12) a velocity asymmetry is therefore less efficient in accelerating the NS than the mass asymmetry in Eqs. (9) and (10). The first-order correction term in s/rs ≪ 1, which we have suppressed when going from Eq. (11) to Eq. (12), reads . It is not necessarily small compared to the leading term , but the corresponding acceleration amplifies the effect of the zeroth-order term, because the displacement by the kick brings the NS closer to the attracting slower parts of the ejecta.

A closer inspection of Eqs. (10) and (12) clarifies the meaning and implications of the non-impulsive nature of the described gravitational tug-boat mechanism for NS acceleration. The instantaneous linear momentum carried by asymmetrically distributed clumps with representative masses and velocities, Pgas ~ 2Δmvs ≈ 2 × 1039 Δm-3vs,5000 g cm s-1, is one to two orders of magnitude smaller than the momentum associated with the estimated possible NS kick, Pns = Mns,1.5vns,8 ≈ 3 × 1041 g cm s-1 with Mns,1.5 = Mns/1.5 M. Efficient NS acceleration can therefore not be explained by considering just the instantaneous momentum of the clumps at any time. Exactly this is the problem connected with purely hydrodynamically mediated kicks where the NS receives its recoil by hydrodynamical forces during the relatively short (hundreds of milliseconds) phase of anisotropic mass accretion and ejection. Correspondingly, Janka & Müller (1994) showed that even with the most extreme assumptions for the explosion asymmetry these more impulsive hydrodynamical kicks are unlikely to lead to recoil velocities of more than a few 100 km s-1. Instead, stronger kicks require a long-lasting accelerating force, which is provided by the gravitational interaction between the NS and the asymmetric ejecta. For this gravitational tug-boat mechanism to work, it is necessary to sustain the momentum of the dense clumps in the ejecta shell over many seconds, despite the gravitational transfer of clump momentum to the NS. This requires a continuous clump acceleration which can happen by hydrodynamical forces within the inhomogeneous ejecta, whereby internal energy of the ejecta is converted into kinetic energy provided the ejecta do not expand ballistically yet. The latter requirement limits the radial scales and time intervals that can contribute to the pulsar kick. The clumps need to be hydrodynamically coupled to the energetic environment they are embedded in. This guarantees that they are continuously accelerated while they transfer momentum and energy to the NS through their long-range gravitational attraction. The energy thus pumped into the motion of the NS can become quite significant,  erg. This, however, is still a small fraction of the  ~1051 erg of a normal SN explosion.

3.5. Progenitor dependence of neutron-star kicks

The data of Table 2 seem to suggest an interesting trend in dependence of the progenitor models: all computed cases for the W15-progenitor yield NS recoil velocities that clearly exceed 200 km s-1, and in all cases the acceleration is still relatively high at the end of the simulation. Also two of the L15 models (the low-resolution run L15-4-lr and L15-5) develop kicks of  ~200 km s-1 and higher with still large acceleration at the end of the simulated evolution, but all models of the B and N progenitors produce clearly slower NSs.

thumbnail Fig. 8

NS kick velocities and acceleration values (indicated by the length of the arrows attached to the symbols) versus explosion energies for all computed models after 1.1–1.4 s of simulated postbounce evolution. Models corresponding to different progenitors are shown by different symbols as defined in the figure.

Open with DEXTER

Both model subsets of strong and weak kicks exhibit a range of explosion energies without any clear correlations between energy and kick velocity (Fig. 8). It is important to note that for each investigated progenitor the number of models is still limited and includes either a wider range of kick velocities for high as well as low explosion energies (W15, L15) or little variation in the kick velocities for all tested energies (B15, N20). We refrain from claiming a lack of high-kick cases for strong explosions (Eexp ≳ 1.2 × 1051 erg) because no W15 simulation was done in this energy range and a trend could be connected with progenitor dependencies. In general, our present results are compatible with conclusions drawn by Scheck et al. (2006) on the basis of a large set of more than 70 simulations in two dimensions and for three different progenitor models with and without rotation. Scheck et al. (2006) observed kick velocities on the low and high sides of the velocity distribution for all investigated progenitors and for less as well as more energetic explosions.

Presently we have no clue which model conditions, i.e., which progenitor or explosion or NS properties, could be favorable or disfavorable for high kick velocities. Neither are W15-models with their typically higher kicks distinct concerning explosion energy or assumed NS contraction, nor is their progenitor structure clearly special in any aspect compared to the other stars. For example, L-models have the highest density and most shallow density profile outside of the iron core and thus the largest mass accretion rates after core bounce with the slowest decay in time. Despite their longest delays until the onset of the explosion, however, they do not stick out concerning their NS kicks.

The determination of clear trends in dependence of progenitor, explosion properties, and remnant properties will therefore require a more extended investigation with a systematic variation in the involved degrees of freedom. This will need many more 3D models than currently available, which is computationally extremely demanding, even with the approximations of the neutrino transport used in our work. At the present time we therefore refrain from drawing premature conclusions or getting entrapped in unfounded speculations.

3.6. Contribution by anisotropic neutrino emission

A directional asymmetry of the neutrino flux radiated by the nascent NS leads to anisotropic energy and momentum loss from the compact remnant. This can provide a kick of the compact remnant, too. The contribution of this neutrino-induced kick to the NS velocity, vns, is equal to the negative of the momentum carried away by the neutrinos divided by the NS mass, (13)This recoil can be expressed in terms of the total neutrino energy loss, Eν, by introducing dimensionless parameters kν,i that describe the asymmetry of the neutrino momentum loss in all three coordinate directions: (14)where c is the speed of light. The vector components kν,i are computed as (15)Here Pij are the components of the neutrino pressure tensor ; Pij defines the flux component j of the neutrino momentum component i. Fν is the neutrino energy flux, dS the surface element on the radiating (neutrino)sphere, (with components nj) the unit vector perpendicular to the surface, and the time and surface integral in the denominator yields Eν. From Eq. (14) an order-of-magnitude estimate of the NS velocity is derived as (16)This yields a recoil velocity of vns ≈ 330 km s-1 for kν = 0.01.

In principle, ultrastrong magnetic fields could lead to large-scale, potentially strong dipolar asymmetries of the neutrino transport out of the hot PNS, e.g. by changing the neutrino-matter interaction opacities, coupling to a neutrino magnetic moment, or affecting neutrino-flavor conversions via the MSW effect (e.g., Bisnovatyi-Kogan 1996; Arras & Lai 1999; Kusenko 2009; Kishimoto 2011). However, little is known about the strength and structure of the magnetic field inside a nascent NS. It could well be chaotic due to convective activity, in which case a large global asymmetry with dominant low-order multipoles is unlikely. Not counting on support by strong dipolar B-field components, a directional bias of the neutrino emission of even only one percent is extremely difficult to obtain.

thumbnail Fig. 9

Neutrino flux asymmetry at  ~150, 300, 600 ms after bounce for model W15-3 as representative case. The 4π-maps show the relative variations ΔF/ ⟨ F ⟩  in the total (i.e., sum for ν and of all three neutrino flavors) neutrino energy flux in all directions. Hot spots (yellow and red) during the postbounce accretion phase (t ≲ 400 ms; top and middle panels) are connected to strong convective downflows in the neutrino-heated postshock layer. The typical scale of the corresponding high-emission areas is 20–30 degrees, and the amplitude of variations ranges around 10–30% during most of the time, though rare, single maxima can reach up to  ~70%. After the accretion has ended (t ~ 400–500 ms), the flux variation pattern is determined by the cell structure of the convection inside the proto-neutron star (i.e., below the neutrinosphere) and the variation amplitute decreases to a few percent (bottom panel). The number in the lower left corner of the plots (dℰ/dt) gives the total neutrino energy loss rate as flux integral over the whole emitting surface.

Open with DEXTER

Asymmetric neutrino emission can be caused by convective overturn in the postshock layer and by convective activity inside the NS, i.e., below the neutrinospheric region. In the postbounce accretion phase, asymmetric neutrino emission is dominated by low-entropy accretion in the neutrino-heating layer between gain radius and shock. Hot spots of high neutrino emission are connected to the strongest accretion funnels. They encompass typical angular scales of 20–30 degrees and are scattered over the whole emitting surface (Fig. 9, top panel). The low-emission regions, which surround the hot spots, trace the distribution of high-entropy bubbles of expanding, neutrino-heated matter and are significantly more extended. The flux variations between low-emission and high-emission regions are usually 10–30%, but in very bright spots can reach up to about 70% for short times (typically of order milliseconds). These local variations are caused by enhanced fluxes of electron neutrinos and antineutrinos, both being abundantly produced in the accretion layer, but the fractional variations are computed relative to the total neutrino flux. The latter includes neutrinos of all three flavors, i.e. also heavy-lepton neutrinos and antineutrinos, whose emission (about 30–50% of the total flux) comes from the core and is not strongly affected by accretion emission and thus more isotropic. After the explosion has been launched, typically a fewer accretion downflows survive and the number of hot spots decreases. The activity is now concentrated in a few regions (Fig. 9, middle panel) instead of showing the previous high-order multipole pattern that encompasses the whole sphere. After the accretion has finally ceased, the surface flux variations are determined by the convectively modified neutrino transport from below the neutrinosphere. The granular pattern of the neutrino emission then reflects the smaller convective cells in the radially less extended convective layer inside the PNS (Fig. 9, lower panel). In this phase the variation amplitude of the total neutrino surface flux is usually around 5%, in peak activity periods up to 10%.

During all three stages, however, the flux variation pattern is extremely time-dependent and non-stationary. The accretion hot spots appear, fade away again, and grow once more at different locations on the surface. Only in the phase of shock acceleration and simultaneous accretion, when the overall structure of the postshock ejecta begins to freeze out and a nearly self-similar expansion sets in, the coherence timescale of the emission pattern increases. Nevertheless, there are still large temporal fluctuations because the accretion flows along the downflows are not continuous but change in strength. In the post-accretion phase the cellular pattern is uniform over the whole surface. The cell sizes show a tendency to increase at later times because the width of the convective layer grows with time. Simultaneously, the amplitude of the flux variations as well as the positions of cooler and hotter areas exhibit rapid temporal variability on timescales of milliseconds to some 10 ms.

As a consequence of the non-stationary, time-variable, and small-scale cellular structure of the surface-flux pattern, in particular during the long-time neutrino emission in the Kelvin-Helmholtz cooling phase of the PNS when most of the neutrino energy is radiated away, we do not expect any sizable NS acceleration due to anisotropic neutrino emission.

thumbnail Fig. 10

Left: entropy distribution in a cross-sectional plane of model L15-1 at 1000 ms after bounce. The plane is chosen such that it contains the final vectors of the NS kick velocity induced by aspherical mass ejection (vns, long white arrow pointing to the right) and of the kick velocity caused by anisotropic neutrino emission (vns, short white arrow pointing towards the 7 o’clock position). The relative lengths of both arrows correspond to the absolute values of the velocities (vns ≈ 160 km s-1 and vns ≈ 5 km s-1), but the short arrow for vns is stretched by a factor of 5 for better visibility. A yardstick in the lower left corner of the figure indicates a spatial length scale of 2000 km. The contribution of anisotropic neutrino emission to the total NS kick velocity is clearly a small effect. Right: time evolution of the neutrino-induced kick velocity vns and of the angle αkν between vns and vns. The final velocities point roughly in opposite directions.

Open with DEXTER

Following Scheck et al. (2006), we calculate the neutrino momentum from the numerical results as (17)where pν =  | Fν | /c2 is the neutrino momentum density, and is the unit vector in radial direction, which in our transport description (“ray-by-ray approximation”) also defines the direction of the neutrino flux (flux components in the angular directions are assumed not to exist). The first term accounts for the total momentum associated with all neutrino energy on the computational grid at a given time t, while the second term is the time integral of the neutrino momentum that has left the computational domain through the outer boundary until this time t.

The numbers for the neutrino-kick velocity vns listed for all models in Table 2 confirm our expectation: the data show that vns makes only a very small, in most cases a negligible, correction to the final NS kick velocity. It is typically much less than 5 km s-1 except in model N20-1-lr, where vns reaches 7 km s-1, which is still only  ~4% of the mass-ejection induced kick velocity vns at 1.3 s after bounce. An entropy slice of model L15-1, which is a better resolved case with a relatively high neutrino-produced kick, is displayed at 1000 ms after core bounce in Fig. 10 together with the time evolution of vns and of the angle αkν between vns(t) and vns(t). The chosen slice contains the final (at t = 1.3 s p.b.) vectors of both velocities. The final angle αkν is significantly larger than 100°, a result that is also found in most of our other 3D models (Table 2).

The velocity vns increases steeply only after about 200 ms p.b. when a sizable accretion and emission asymmetry begins to develop. Before this phase the small angular variations in the neutrino emission do not allow for the growth of any noticeable value of vns. Towards the onset of the explosion (at  ~400 ms; Table 2) a nearly stationary pattern of accretion downflows and high-entropy bubbles emerges in the postshock layer and persists during the subsequent expansion of the supernova shock and neutrino-heated ejecta. The steep rise of vns between 200 ms and 600 ms after bounce in model L15-1 is associated with a massive and extended, long-lasting downflow funnel that carries fresh gas swept up by the accelerating shock from the immediate postshock region towards the NS surface. This leads to a neutrino emission maximum in the right hemisphere (roughly at the 2:30 o’clock position) and thus to a neutrino-induced kick pointing in the opposite direction (short white arrow for vns in Fig. 10). The relics of the massive downflow funnel are still visible to the right of the center in the left image of Fig. 10. At about 600 ms p.b. the accretion ends and the neutrino emission pattern adopts the high-multipole cellular structure characteristic of the PNS convection during the subsequent Kelvin-Helmholtz cooling phase (cf. above discussion and bottom panel in Fig. 9). At this stage the neutrino-kick velocity has nearly reached its terminal value. It hardly changes later on because short-wavelength variations and short-timescale fluctuations diminish the integrated neutrino emission asymmetry over longer periods of time. Nevertheless, αkν continues to exhibit some drift even between 600 ms and  ~1000 ms, because the gravitational pull by the anisotropic ejecta still causes an evolution of the associated kick direction. Only at t ≳ 1000 ms vns has also found its final orientation about 135° away from vns (Fig. 10), while the gravitational acceleration of the NS continues on a significant level even beyond 1.3 s after bounce (Table 2). The neutrino kick typically damps the ejecta-induced kick: large angles αkν are a natural consequence because long-lasting downflows determine the accelerating gravitational tug on the one hand and produce the dominant neutrino emission anisotropy on the other. While the NS kick associated with the asymmetric ejecta distribution is directed towards the downflow, the enhanced neutrino emission on the downflow side carries away momentum that leads to a NS recoil in the opposite direction.

4. Neutron star spins

Mechanisms that impart a natal kick to the NS can also lead to a NS spin even in the case of nonrotating progenitor stars. This is expected theoretically (Burrows et al. 1995; Spruit & Phinney 1998; Fryer & Young 2007) and suggested by observations (e.g., Farr et al. 2011). NS spin-up may, for example, be a stochastic process that is connected to random, nonradial impacts of the fast accretion downflows, by which low-entropy gas is channelled through the gain region onto the PNS. Accretion downflows, which hit the NS asymmetrically and not exactly head-on, exert torques and thus spin up the NS. A mass Δm colliding with the NS with impact velocity and impact parameter d transfers an angular momentum of (18)and for rigid rotation leads to a spin period of (19)where Δm is normalized by 10-3M, d by 30 km, vi is given in units of 1010 cm s-1, and Ins ~ 2  ×  1045 g cm2 is the NS moment of inertia.

Equations (18) and (19) give rough order of magnitude estimates and show that for moderate values of the nonaxisymmetrically accreted mass (Δm ~ 10-3M) the NS may receive an angular momentum of some 1046 g cm2 s-1, corresponding to a specific angular momentum of order 1013 cm2 s-1 and spin periods of some 100 ms. However, neither the number of accretion downflows nor their impact parameters can be computed analytically. Quantitative answers require 3D hydrodynamical simulations of the postbounce accretion phase.

Alternatively, it has been proposed that NS spin-up may not result from a stochastic process connected to random, nonradial impacts of accretion downflows. Instead, SASI spiral modes (m = 1, 2 spherical harmonics components) might develop in the postshock accretion flow, by which a radial redistribution and separation of angular momentum between the compact remnant and the ejecta could be established: While NS spin is accumulated by the accretion of rotating gas, the ejecta carry away counter-rotating matter with an angular momentum of the same magnitude but with opposite direction (Blondin & Mezzacappa 2007; Blondin & Shaw 2007; Iwakami et al. 2009; Fernández 2010). This phenomenon was seen in idealized numerical models of collapsing stellar cores, making use of an ideal-gas EoS and ignoring neutrino heating (except by Iwakami et al. 2009). It was also observed in laboratory experiments of a shallow water analogue of the SASI in two dimensions (Foglizzo et al. 2012). However, the timescale for SASI spiral modes to develop in the core without initial rotation was found to be close to one second, which is longer than the timescale for the onset of the explosion and for the saturation of the NS spins in all of our SN simulations. Correspondingly, we do not observe any clear signatures of SASI spiral modes in our models. In contrast, m = 1, 2 modes can dominate the asphericity of the shock surface, Rs(θ,φ), and of the postshock mass distribution, (Eq. (3)), in models where the NS acquires only an unspectacular angular momentum. In all our models, however, the normalized amplitudes of the m ≠ 0 spherical harmonics components of the shock deformation never exceed a few percent (Fig. 5). This is much smaller than the order unity shock displacements obtained by Blondin & Mezzacappa (2007); Blondin & Shaw (2007) and Fernández (2010). It is also possible that neutrino heating and neutrino-driven convection destroy the coherence of the spiral modes and thus impede their growth, or it might enhance mode coupling, thus amplifying energy redistribition between different modes. Moreover, even when shells with strong differential rotation in opposing directions arise in the turbulent region between PNS and shock, the question still remains how much of the induced angular momentum is finally accreted. Rantsiou et al. (2011) suspect that most of it is likely to be ejected in the explosion.

thumbnail Fig. 11

Time evolution of the NS kick velocity (upper left panel) compared to the absolute value of the NS angular momentum (lower left panel), the polar and azimuthal angles (θspin, φspin) of the NS spin direction (right top and middle panels), and the relative angle αsk between spin and kick directions (right bottom panel) for five selected models. Two of these models represent high-spin cases (N20-2, N20-3), two are low-spin cases (W15-6, W15-7), and one is an intermediate case (L15-3) from our model set. The build-up of the angular momentum occurs in mainly a single impulsive episode during the convectively perturbed accretion phase of the NS. Together with its angles the angular momentum reaches its terminal value long before the NS acceleration has ended. The clustering of all five plotted models around a final θspin ≈ 90° is a mere coincidence, and the whole set of computed models exhibits the expected random distribution in polar-angle space.

Open with DEXTER

4.1. Numerical results

Because of violent convective activity, the matter in the neutrino-heated layer develops nonradial motion in different directions. These nonradial flows can be connected with sizable amounts of angular momentum. Integral values can indeed reach up to several 1046 g cm2 s-1 in different regions of the layer between PNS surface and accretion shock, compatible with the order-of-magnitude estimate of Eq. (18).

Since our simulations start from nonrotating collapsing stellar cores, any angular momentum acquired by the nascent NS must be balanced by angular momentum in the opposite direction in the gas outside of the NS. Different from the linear momentum of the compact remnant, which is calculated as the negative value of the gas momentum integrated over the whole computational grid (see Eq. (1)), numerical tests showed that it is advisable to constrain the volume integration for the angular momentum to the immediate vicinity of the forming NS. This helps minimizing numerical errors connected to the imperfect conservation of angular momentum of gas that expands outwards over radial distances extending over many orders of magnitude on a spherical grid with varied resolution. As in Wongwathanarat et al. (2010b) we therefore determine an estimate of the angular momentum Jns transferred to the NS by hydrodynamic forces associated with accreted or ejected gas from the following equation: (20)where j is the specific angular momentum, Ω denotes the solid angle, and is the area of a sphere of chosen radius ro. The first term in the brackets on the rhs of Eq. (20) accounts for the angular momentum contained by the spherical shell bounded by the NS “surface” at Rns (defined as the radius where the density is 1011 g cm-3) on the one side and bounded by ro on the other. The second term gives the angular momentum that is carried by mass leaving the shell through surface Ao until time t. The introduction of the surface integral is possible in the context of NS spin-up, because angular momentum can be transferred to the central object only when the latter exchanges mass with its environment. Long-range gravitational forces cannot exert a torque on the compact remnant, because the NS (in a very good approximation) is described in our simulations as a gravitating point mass plus an essentially spherical near-surface layer of matter surrounding the inner grid boundary. The numerical evaluation of Eq. (20) yields results that do not depend on the exact location of ro between 500 km and 1000 km.

In Table 2 the results for the angular momentum Jns, normalized by 1046 g cm2 s-1, the relative angle between spin and kick directions, αsk, and the estimated final NS spin period, Tspin, are listed for all our model runs at the end of most of the simulations (1.1–1.4 s after bounce). We compute Tspin from the given numbers for Jns and the NS moment of inertia, Ins, as Tspin = 2πIns/ | Jns | , considering the NS as a rigid rotator and making the rough approximation that Ins is given by the Newtonian expression for a homogeneous sphere: . The (baryonic) NS mass Mns is also listed in Table 2 and we adopt a final NS radius of Rns = 12 km. It is important to note that we assume that the angular momentum of the NS is conserved when the remnant contracts to its final radius during the neutrino-cooling evolution. This ignores a possible  ~40% reduction associated with neutrino losses (see Janka 2004) and with mass ejection in the neutrino-driven wind.

The angular momentum values listed in Table 2 span a range from  ~0.1 × 1046 g cm2 s-1 to more than 7 × 1046 g cm2 s-1, while the corresponding spin periods vary by over a factor of 60, ranging from close to 0.1 s to nearly 8 s. Higher angular momenta usually lead to the smaller Tspin; this order is reversed by differences of the NS masses only in a few cases.

While there is no obvious systematic differences of the NS angular momenta between the W15 and L15 models, although the L models clearly have the longest delays to explosion, three of the four members of the N15 series stick out having the highest Jns values of the whole model set. In addition, two of the three B models are among those acquiring the lowest NS angular momentum. These tendencies, cautiously interpreted in view of the limited number of models in each series, seem to be linked to the convectively stirred ejecta mass in the neutrino-heating region. Assuming similar values for the specific angular momentum j in all cases, the integration of Eq. (20) depends crucially on the mass involved. In particular model B15 is special by its steep density decline outside the iron core in addition to the short evolution timescale to explosion, both of which tend to give a small angular momentum transfer to the accreting remnant.

The values of the final angle between spin and kick vectors, αsk, in Table 2 do not reveal any clear correlations beween spin and kick directions, in particular no strong tendency towards alignment or counter-alignment. The distribution is nearly symmetric for αsk < 90° and αsk > 90° (9:11 cases) and essentially flat when sorted into 60° intervals (7:6:7 cases), although a statistically equal distribution in all directions would require twice as many events in the equatorial belt than in each polar wedge. Our model sample, however, is still too small to draw firm conclusions in a statistical sense.

4.2. Neutron star spin-up mechanism

It is important to note that kicks and spins are attained by the NS through different mechanisms in our simulations, although both phenomena are connected to the development of large-scale nonradial flows in the collapsing stellar core, which lead to momentum and angular momentum exchange between the accreting compact remnant and the asymmetrically expanding ejecta gas. In order to understand the differences, we present in Fig. 11 the time evolution of Jns compared to vns for a selection of representative models, namely for the two high-resolution N15-runs with the largest values of Jns, for the two W15-cases with the lowest Jns, and for model L15-3 as an intermediate case. The left two panels show that the angular momentum increases much faster and reaches a value close to its asymptotic one much earlier than the recoil velocity does. Different from the NS kick most of the NS spin is attained in an impulsive event within a time interval of 100–200 ms relatively early after core bounce. No significant Jns evolution takes place beyond t ≳ 500 ms, i.e. after the onset of the explosion (measured by texp). In contrast, the NS velocity is very small at the time when the PNS receives its angular momentum, but the kick can grow monotonically with considerable rate for several seconds.

The reason for this difference is the fact that angular momentum can be transferred to the NS only by nonspherical mass flows and thus during a few hundred milliseconds of accretion. This takes place before the onset of the explosion and can go on afterwards only as long as accretion still continues simultaneously with the outward expansion of the SN shock5. In contrast, the recoil of the newly formed NS is mainly established by the long-range anisotropic gravitational attraction between the compact remnant and the asymmetric ejecta as described in detail in Sect. 3. In Fig. 11 this difference manifests itself in the long-lasting, monotonic rise of vns(t) during the whole displayed time and even much longer (cf. Fig. 4), whereas most of the angular momentum is accumulated in the NS during a period of massive accretion shortly after the convective and SASI instabilities have grown to their full strength. When the accretion to the NS is diminished by the outward acceleration of the gas behind the explosion shock, sporadic downdrafts reaching the NS can still lead to modest variations in Jns for some 100 ms. The duration of this phase depends on the progenitor structure and the detailed dynamics of the explosion, but even for relatively low explosion energies it does not last any longer than  ~500 ms after the revival of the SN shock.

The differences between the kick and spin mechanisms are reflected not only by the evolution of the NS kick and spin magnitudes but also by the evolution of the kick and spin directions. During the initial steep growth of Jns the spin direction (θspin, φspin) mostly settles to its final direction, in particular in the two high-spin N20-models. Long-time variations in the NS angular momentum then correlate with only smaller changes in the spin direction, which signal the impact of new accretion downdrafts. Such changes become less frequent at later times and finally die away when the accretion ends. In the L15-models, which have the highest densities in the shells surrounding the iron core and therefore explode with the longest delays (a good example is model L15-3, where the shock is revived as late as  ~500 ms after bounce), the period of large Jns variations is more extended. Slow fluctuations of the spin-vector angles and a slight drift of Jns occur even until more than one second post bounce because of the angular momentum transported away with the more massive neutrino-driven wind in the L models, in contrast to the W and N models. The long-time variation in the spin differs from that of the kick, whose direction freezes out as soon as the bulk of the mass carrying the ejecta asymmetries expands essentially self-similarly, although vns continues to increase. The bottom right panel of Fig. 11 provides the time evolution of the relative angle between kick and spin vectors. The high-frequency variations in αsk are caused by rapid changes in the kick direction (θkick, φkick) until about 200 ms after the onset of the explosion, while a slow drift of αsk on longer timescales (up to  ~1 s in the most extreme case) is connected to an ongoing reorientation of the spin direction.

Interestingly, the highest NS kicks and spins are not obtained for the same models or progenitor. While the remnants of the N20 progenitor exhibit the tendency to develop the highest spins, those of the W15 series are found to end up with the largest kicks. This can be considered as another manifestation that spin and kick mechanisms are linked to different aspects of anisotropic mass flows in the SN core, namely the former to asymmetries of the accretion flow and the latter mainly to asphericities of the ejecta distribution.

It is unclear whether NS spin periods that are more extreme than those found in some of our N20 models, in particular periods below 100 ms, can be obtained in nonrotating stars. It is possible that faster NS rotation requires the progenitor core to rotate (e.g., Heger et al. 2005; Ott et al. 2006). In a rotating environment, however, the growth of SASI spiral modes will be severely altered and potentially fostered (Blondin & Mezzacappa 2007; Yamasaki & Foglizzo 2008; Iwakami et al. 2009), and three-dimensional simulations are therefore indispensable to make predictions for the connection between progenitor and NS rotation.

5. Heavy-element production and explosion asymmetries

Fryer & Kusenko (2006) pointed out that ejecta asymmetries could be used to distinguish between different theoretical suggestions for the kick mechanism. Ejecta-driven kicks like the gravitational tug-boat mechanism discussed here must be expected to be associated with stronger ejecta motion, i.e. with higher ejecta momentum, in the direction opposite to the kick. However, diagnosing momentum asymmetries of the ejecta may be difficult except in cases where the center of mass of the gaseous SN remnant is clearly displaced relative to the NS position or where the gas cloud exhibits extreme deformation without gradients of the environmental conditions being responsible for the asymmetric ejecta expansion. The SN remnants Cassiopeia A (Isensee et al. 2010; DeLaney et al. 2010; Rest et al. 2011; Hwang & Laming 2012) and Puppis A (Petre et al. 1996; Winkler & Petre 2007; Katsuda et al. 2008, 2010; Becker et al. 2012) may be such lucky cases, because in the former all the SN ejecta in the remnant, whose mass is dominated by oxygen, while in the latter fast-moving, oxygen-rich optical filaments and knots seem to have been recoiled opposite to the direction of the high-velocity compact stellar remnant. In general, however, a major problem for observing momenta asymmetries linked to an asymmetric beginning of the explosion may result from the fact that the outgoing SN shock transfers the bulk of the explosion energy and ejecta momentum to the He and H shells, which dominate the ejecta mass in events other than stripped SNe. An initially asymmetric explosion thus becomes much more spherical when the SN shock sweeps up massive stellar He and H layers. In this context it is important to note that the NS kinetic energy and the NS momentum are very small compared to the kinetic energy of the SN explosion, Ek,SN, and the corresponding measurable momentum, , associated with the radial motion of the ejecta mass Mej. To determine the (relatively small) linear momentum of the ejecta, which balances the momentum of the NS, it is necessary to sum up all components of the momentum vectors of the whole ejecta mass. Because of projection effects this may be difficult and error-prone even if all of the ejecta can be observed.

For this reason we propose here to look for correlations between NS kicks and asymmetries in the distribution of heavy elements. These are explosively produced in the innermost, shock-heated (and neutrino-heated) SN matter during the first seconds of the blast, and their globally asymmetric distribution will not be destroyed lateron. In particular iron-group nuclei, most of which are radioactive 56Ni decaying to stable 56Fe, are good tracers of such early explosion asymmetries. This is connected to the fact that the initial Fe-group material in the progenitor core is completely photo-disintegrated during collapse and then buried in the forming NS. All expelled Fe-group nuclei are freshly assembled in the neutrino-heated SN ejecta and in the silicon and oxygen layers that are heated to sufficiently high temperatures by the accelerating SN shock.

thumbnail Fig. 12

56Ni distribution for model W15-2 at 365 ms and 515 ms after bounce as produced by explosive burning in shock-heated ejecta. The integrated nickel mass per solid angle is plotted in the cross-sectional plane that contains the NS kick (indicated by the black arrows pointing to the right) and spin vectors. Red colored regions show that more nickel is formed in the direction opposite to the NS kick direction. The size of the displayed volume can be estimated from the yardsticks given in the lower left corner of each image.

Open with DEXTER

Figure 12 visualizes the corresponding directional asymmetry of the 56Ni production in the plane containing the NS velocity and spin vectors in the high-kick model W15-2. The images show the nickel mass per unit solid angle at post-bounce times of 365 ms and 515 ms (the explosion in this model sets in about 250 ms after bounce; Table 2), which is obtained by integrating the nickel-mass density over radius within each zone: (21)where XNi(r) is the mass fraction of nickel and ρ(r) the mass density at location r. The color coding clearly shows more nickel production (intense red) in the direction opposite to the kick vector.

At around 450–500 ms after bounce the nickel fusion in the shock-heated layers comes to an end. While most of the nickel expands in a thin and dense shell behind the shock, some of it is channeled into massive clump-like “pockets” between the more dilute high-entropy bubbles of rising and outward pushing neutrino-heated matter. Also these bubbles may contain a fair amount of nickel from the α-rich freeze-out that occurs during the expansion cooling of the bubbles. However, because these regions have a slight neutron excess in our models, the nuclear flow is assumed to lead to the production of the “tracer” species included in our network as described in Sect. 2.3. Since the actual composition of this tracer material strongly depends on the inaccurately determined value of Ye set by νe and absorption and emission in the neutrino-processed ejecta, the amount of 56Ni in this material is uncertain and we have not depicted it in Fig. 12. The tracer material, however, enhances the hemispheric asymmetry of the nickel distribution. This is visible from the data listed in Table 3, which provides abundance yields of the nuclear species in our burning network after all nucleosynthesis processes have ended (depending on the progenitor and explosion energy, this is the case typically 100–300 s after the onset of the explosion). Higher values of 56Ni in the southern hemisphere correlate with a similarly large excess of tracer material on the same side. We define the northern hemisphere as the half-sphere around the NS kick direction, which implies that more 56Ni and tracer material are ejected opposite to the NS kick direction. This effect is especially pronounced (with hemispheric contrast up to a factor 3.5) in high-kick models like W15-1 and W15-2. The last column of Table 3 provides our determination of the total 56Ni yields of the listed models. The uncertainty range is defined by the assumption that in the lower limit the tracer material contains no 56Ni and in the upper limit all tracer material is 56Ni.

thumbnail Fig. 13

Shock deformation for the high-kick models W15-1 and W15-2 and for the moderate-kick models L15-1 and L15-2 at 1.3 s or 1.4 s after bounce, as given in the upper right corner of each map. These times are close to the evolution phases when the nickel production in the shock-heated ejecta takes place. The color scale is linear with black indicating minima and red maxima. The differences between maxima and minima relative to the average shock radii are 30%, 35%, 34%, and 29%, respectively. Red crosses indicate the NS kick directions while white plus signs mark the directions opposite to the kicks. The white lines separate the corresponding two hemispheres. The northern hemisphere is defined as centered around the NS kick vector.

Open with DEXTER

Table 3

Hemispheric ejecta yields of helium and metals, and total nickel mass for the high-kick models W15-1 and W15-2, and the moderate-kick models L15-1 and L15-2.

thumbnail Fig. 14

Integrated nickel masses per solid angle in all directions displayed in 4π maps for models W15-1, W15-2, L15-1, and L15-2 at postbounce times well after all nucleosynthesis processes have seized in our simulations. Red crosses indicate the NS kick direction while white plus signs mark the direction opposite to the kick. The corresponding hemispheres are separated by a white equator line. Note the differences between the high-kick W15 models and the moderate-kick L15 models. The former show a clear enhancement of the nickel production in the hemisphere opposite to the kick, whereas the latter eject nickel significantly more isotropically.

Open with DEXTER

thumbnail Fig. 15

Volumetric three-dimensional visualization of the nickel distribution for models W15-1, W15-2, L15-1, and L15-2 at the postbounce time given in the top left corner of each panel along with the model name. The semi-transparent isosurfaces correspond to a chosen value of the nickel mass per grid cell of 3 × 1026 g and are displayed at a stage well after all nucleosynthesis processes have seized in our simulations. The orange vectors represent the NS kick directions and are scaled by the corresponding NS velocities (96 km s-1 for the shortest arrow and 575 km s-1 for the longest one). The high-kick models W15-1 and W15-2 in the upper two panels exhibit a clear asymmetry with much more nickel being ejected in the hemisphere opposite to the kick direction. In contrast, the moderate-kick models L15-1 and L15-2 in the lower two panels exhibit a more isotropic distribution of the nickel, in particular no obvious hemispheric asymmetry between kick and anti-kick directions. While the radial distribution of the nickel may be strongly affected and changed by subsequent mixing instabilities that develop after the outgoing shock has passed the composition-shell interfaces of the progenitor star, the hemispheric differences in the nickel ejection will not be destroyed during the later supernova explosion.

Open with DEXTER

The comparison of the high-kick models W15-1 and W15-2, both of which have NS recoil velocities in excess of 500 km s-1 (Table 2), with the moderate-kick model L15-1 (vns ~ 250 km s-1) and model L15-2, whose NS receives a recoil velocity of only  ~100 km s-1, is enlightening. The data in Table 3 demonstrate that the hemispheric differences of nickel and tracer ejection are largest for W15-1 and W15-2, still considerable for L15-1, and nearly disappear for L15-2. Similarly large hemispheric differences as for 56Ni and our tracer species show up in the masses channeled by the network into 28Si, 40Ca, and 44Ti. All of these nuclear species6 emerge from the innermost ejecta regions, which carry the explosion asymmetries and the mass-distribution inhomogeneities that are responsible for the gravitational acceleration of the newly formed NS. While some degree of hemispheric differences correlating with the NS kick magnitude is still visible for 24Mg, the effect in the case of 16O is very small. Although oxygen is more abundantly synthesized behind the stronger shock in the southern hemisphere, an asymmetry of the distribution of this element is largely obscured by the dominance of the spherically symmetric ejecta mass that originates from the unperturbed oxygen shell of the progenitor star. The other lighter species (20Ne, 12C, and 4He) do not exhibit the southern excess. In contrast, their yields are slightly reduced in the southern hemisphere – most pronouncedly in the case of 20Ne – because of the explosive production of the higher-mass nuclei at the expense of these lighter abundances.

The dependence of the explosion strength on direction is visualized in Fig. 13 in terms of the variations in the shock radius and thus shock velocity at around the time when the nickel nucleosynthesis has just taken place. Figure 14 shows the corresponding radially integrated nickel masses per solid angle, (22)where the radial integration is performed between the inner and outer grid boundaries, Rib and Rob, respectively. In Fig. 15 we provide volumetric impressions of the nickel distribution in the ejecta by plotting isosurfaces of a suitably chosen constant value of the nickel mass per grid cell7. The coordinate system in these images has been rotated such that the NS kick vector points to the 12 o’clock position.

In Fig. 13 one can see clear differences of the shock deformation between low-kick and high-kick models. In the high-kick cases of W15-1 and W15-2, the northern side, which is defined as the hemisphere around the NS kick direction (marked by the red cross), contains a huge, deep trough (blue and black) whose minimum roughly coincides with the kick direction. In contrast, the southern hemisphere is filled by a wide bump with several local maxima. Consequently, there is a clear dipolar shock deformation along the kick axis (cf. also Fig. 5). The moderate-kick model L15-1 still exhibits a tendency for the same asymmetry pattern with a shock minimum close to the red cross and a maximum near the white plus sign in the opposite direction. Now, however, shallow minima and maxima also appear in the southern and northern hemisphere, respectively, indicating the growth of a quadrupolar deformation component. The low-kick model L15-2 violates the north-south imbalance. Both shock maxima and minima cross the equatorial line and are obviously not correlated with the locations of the cross and plus signs. Instead of the dipolar shock deformation of models W15-1 and W15-2 we now see a quadrupolar pattern, consistent with the left lower panel of Fig. 5. We note that the relative variation in the shock radius in all four displayed cases is similar. The shock deformation in high-kick and low-kick models therefore does not need to differ in the variation amplitude, but instead the type of dominant spherical harmonics mode can be the more important factor.

In Fig. 14 the highest concentrations of ejected nickel are found exactly at the locations of the shock maxima in Fig. 13, while big “holes” of the nickel distribution appear in the regions of shock minima. Again, the angular pattern of the 56Ni distribution of the low-kick model L15-2 is clearly dominated by higher spherical harmonics modes than in the high-kick models W15-1 and W15-2. Model L15-1 with a moderate kick defines an intermediate case. The directional asymmetries depicted in Fig. 14 are complemented by the volumetric information of Fig. 15. The black areas seen in Fig. 14 correspond to directions where holes exist in the spatial nickel distribution, while the green-yellow-red maxima in Fig. 14 appear as big, clump-like, coherent structures of nickel in Fig. 15. The hemispheric nickel imbalance in models W15-1 and W15-2 is obvious, whereas L15-1 shows a bipolar structure that is misaligned with the NS kick vector. Different from the three others, L15-2 exhibits a more isotropic nickel distribution and considerably more power in structures on smaller angular scales.

In summary, Figs. 1215 confirm that the production of nickel and of adjacent products of explosive burning (cf. Table 3) are excellent tracers of explosion asymmetries. An observational determination of major anisotropies of the motion of heavy elements between silicon and nickel in SN remnants and an expansion of these metals mainly opposite to the direction of the NS motion would lend strong support to the NS acceleration by the described gravitational tug-boat mechanism.

6. Summary and conclusions

We have performed a set of 20 SN simulations from shortly after core bounce to more than one second after the onset of the explosion for four different 15 M and 20 M nonrotating progenitor stars, making use of an axis-free Yin-Yang grid in three dimensions. The grid helped us to avoid numerical artifacts as well as severe time-step constraints near the polar axis of a usual spherical grid. A further gain in performance was achieved by excising the high-density core of the PNS (at densities significantly above the neutrinosphere) and replacing it by an inner boundary condition, where the contraction of the PNS core was imposed and the neutrino luminosities were chosen such that neutrino-energy deposition behind the shock led to SN explosions of defined energy. The neutrino transport on the grid was treated by a computationally efficient, analytic characteristics integrator, in which a grey approximation by was used averaging over Fermi-Dirac spectra and flux variations with latitude and longitude were represented by employing a ray-by-ray description (see Scheck et al. 2006).

A subset of our explosion models was evolved beyond three seconds postbounce to follow the approach of the NS recoil velocity to the saturation level. Some of these simulations were carried on for several 100 s to determine the explosive nucleosynthesis of heavy elements and the spatial distribution of the burning products.

Our results demonstrate the viability of the gravitational tug-boat mechanism for NS acceleration in three dimensions. In this scenario the newly formed NS is accelerated over a timescale of many seconds mainly by the anisotropic gravitational attraction of the asymmetrically expanding neutrino-heated ejecta gas, in which hydrodynamic instabilities had created large-scale asphericities prior to the onset of the explosion. Before the accretion flows to the NS are quenched by the accelerating explosion, hydrodynamic forces on the NS play a nonnegligible role, too. We emphasize, however, that anisotropies of the neutrino emission, which our ray-by-ray transport approximation tends to overestimate, contribute to the NS kick velocity only on a minor level of percents or less. The reason for this small effect is simple: anisotropies of the neutrino emission exhibit short-timescale intermittency and are characterized by a higher-order multipole structure in angular space. As a consequence, the neutrino-induced momentum transfer is diminished by statistical averaging and is therefore unable to cause any significant net kick in a certain direction. Moreover, only a small fraction of the total neutrino energy loss is radiated anisotropically during the simulation time. Instead, the far majority of the escaping neutrinos contribute to the isotropic background flow that diffuses out of the essentially spherical, hot accretion layer around the nascent NS.

NS acceleration by the gravitational tug-boat mechanism was first discussed on the basis of axisymmetric (2D) simulations by Scheck et al. (2004, 2006) and was confirmed in 3D with a small set of explosion models by Wongwathanarat et al. (2010b). Moreover, support for this mechanism was recently also provided by the 2D simulations of Nordhaus et al. (2010, 2012), in which the NS was allowed to move self-consistently out of the grid center.

While the artificial constraint to axisymmetry enforces a collimation of large-scale flows parallel to the polar grid axis and thus tends to favor a pronounced, dipolar deformation of the explosion, the 3D asphericities appear less extreme and seemingly less promising for high NS kicks. Nevertheless, the inhomogeneities and anisotropies in the massive postshock shell are sizable also in 3D. In fact, low spherical harmonics modes of dipolar and quadrupolar character can dominate the asymmetry of the shock and postshock ejecta also in 3D, although with somewhat smaller amplitude than in extreme 2D situations (Fig. 5). Correspondingly, also in our set of 20 three-dimensional models several cases were found that produced NS kick velocities well above of 500 km s-1, and our current record holder developed a NS recoil of more than 700 km s-1 after 3.3 s (Table 2 and Fig. 4). This is compatible with the 2D results of Scheck et al. (2006), who obtained velocities up to 800 km s-1 after one second and estimated the final NS speed in one out of more than 70 models to be nearly 1000 km s-1 and in another one to be more than 1200 km s-1. Simple analytic estimates of the long-time influence of the gravitational forces by the expanding SN ejecta confirm that dense clumps of some 10-3M (accounting for mass distribution asymmetries of around 10% in the postshock ejecta shell) could well accelerate the NS to velocities of 1000–2000 km s-1 as measured for the fastest young pulsars (cf. Eq. (10)). With an ejecta-shell mass of a few 10-2M, a similar accelerating effect on the NS can be achieved by hemispheric asymmetries of the expansion velocity of several 10% (cf. Eq. (12)), which is in the ballpark of shock expansion differences seen in some of our models (Fig. 13).

The long-time acceleration of the NS is therefore connected to slow-moving, massive ejecta “clumps”, whose expansion lags behind the rest of the ejected material of the SN core. For this reason the NS gains momentum opposite to the direction of fastest expansion of the explosion shock and postshock shell. Correspondingly, the explosive nucleosynthesis of heavy elements from 28Si to the iron group, especially also of radioactive 56Ni, occurs primarily in the hemisphere pointing away from the NS velocity vector. This effect clearly correlates in strength with the size of the kick. We therefore propose that an observational confirmation of opposite momentum directions of the NS and of the bulk of elements heavier than silicon in supernova remnants would lend support to the gravitational kick mechanism discussed in this paper. Such an analysis for the Cassiopeia A (Cas A) SN remnant with its central point source could be enlightening.

Nonradial (off-center) accretion downflows in the postshock layer can also lead to a spin-up of the nascent NS in nonrotating progenitor stars. This requires angular momentum separation in the postshock flow during the postbounce accretion phase of the PNS. The corresponding angular momentum can accumulate to values of several 1046 g cm2 s-1 in our simulations, in agreement with an order-of-magnitude estimate on grounds of simple considerations (Eq. (18)). Assuming angular momentum conservation during the later evolution, we derive spin periods between  ~100 ms and several seconds for the NS in our models. We note, however, that angular momentum conservation is difficult to assure in the numerical models, and thus these numbers should be considered only as a rough range of possibilities.

While the NS acceleration is a long-time (continuing over seconds) phenomenon mostly achieved by the gravitational pull of the postshock ejecta shell, the spin-up of the nascent NS is connected to anisotropic mass flows that transfer angular momentum to the accretor. The relevant postbounce accretion phase of the PNS may last for several 100 ms beyond the onset of the SN blast. This limits the time interval for the spin-up of the NS to a much shorter period than the NS acceleration. Kick and spin of the NS in this scenario are therefore not obtained by the same process: The former is linked to mass ejection asymmetries, the latter to anisotropic, nonradial accretion flows; they are mediated by different forces – gravitational versus hydrodynamical –, and they also build up on different timescales. It is therefore not astonishing that we do not find any correlation of spin and kick directions in our explosion models of nonrotating progenitor stars, and a possible spin-kick alignment suggested by observations (e.g., Johnston et al. 2005, 2007; Kaplan et al. 2008; also see Wang et al. 2006, 2007) remains unexplained.

In this context, the question of the consequences of (even a modest amount of) angular momentum in the progenitor core is interesting and deserves future exploration. In the described kick scenario, however, rotational averaging of the kick, which could cause an alignment of pulsar velocity and spin as suggested by Lai et al. (2001) and Wang et al. (2007), is unlikely to be efficient despite the long timescale of the NS acceleration. Rotational motions in the expelled gas, whose clumpy asymmetries accelerate the NS by gravitational forces, are diminished during expansion, because under the constraint of angular momentum conservation the angular velocity decays with growing radius like r-2. It is therefore not possible to destroy the randomness of the kick direction just by the rotation of the accelerating ejecta asymmetries. An asphericity that tends to accelerate the NS along the spin direction would have to be imprinted in the ejecta gas already at the onset of the explosion. It is unclear whether core rotation can produce such an asymmetry in the case of reasonable and realistic assumptions about the spin rate of the progenitor core prior to collapse. With predicted typical Fe-core rotation periods of (Heger et al. 2005), which lead to NS spin periods of , the neutrino-heated postshock layer has rotation periods of about one second. This is much too slow to have any dynamical influence. On the other hand, if sufficiently fast rotation allowed for a rapid growth of SASI spiral modes or triaxial dynamical or secular instabilities, kicks perpendicular to the NS spin appear more likely than the alignment of both. In our simulations we do not find any evidence of the large-amplitude SASI spiral modes that were observed in idealized setups (without neutrino effects) by (Blondin & Mezzacappa 2007) and Fernández (2010) and experimentally by Foglizzo et al. (2012). The growth of such large-scale, coherent patterns in the postshock accretion flow might be inhibited by the presence of neutrino heating, or their growth timescale might be longer than the evolution period until the onset of the explosion in our models. Numerical studies (Blondin & Mezzacappa 2007; Iwakami et al. 2009) and analytical treatment (Yamasaki & Foglizzo 2008), however, predict faster growth of SASI spiral modes in rotating environments.

Our still limited set of simulations does not allow us to draw any firm conclusions on progenitor or explosion dependences of the spin and kick magnitudes. One of the employed 15 M stars (the W15 progenitor) seems to exhibit a tendency to favor higher recoil velocities on average than the other three investigated progenitors, but the exact reason for this trend is not determined.

Podsiadlowski et al. (2004) hypothesized that the collapse and explosion of low-mass SN progenitors with ONeMg cores leads to NS with low-velocity kicks. They argue that in the case of ONeMg cores in contrast to more massive Fe-core progenitors the explosion sets in more quickly and does not provide the time for the growth of large-amplitude asymmetries of low-order multipole type in the neutrino-heated accretion layer between NS and stalled shock. Indeed, the bounce shock accelerates quickly in the steep density profile at the edge of collapsing ONeMg cores and the explosion sets in on a timescale of only  ~100 ms after core bounce (Kitaura et al. 2006; Janka et al. 2008). This leaves time for only one convective turnover. Therefore, the angular size of the convective plumes remains rather small. Because of the extremely fast, spherical expansion of the shock, the convective structures are also rather symmetric in all directions (cf. Fig. 1 in Wanajo et al. 2011). These facts together with the extremely small mass swept up from the dilute environment of the ONeMg core by the outgoing shock disfavor high NS kicks by the gravitational tug-boat mechanism8. In our scenario of NS acceleration the differences of low-kick and high-kick cases are connected to differences in the explosion dynamics, manifesting themselves in a different growth of hydrodynamic (convective and SASI) instabilities during the onset of the explosion. The growth of these instabilities can be seeded by small random perturbations (which should be present in any convectively perturbed stellar environment and are imposed by us artificially on spherical pre- or post-collapse progenitor models) and leads to anisotropic mass ejection. Large-scale asymmetries that might develop in the progenitor core during convective shell-burning phases prior to collapse as advocated by Burrows & Hayes (1996); van den Heuvel (2010), and Arnett & Meakin (2011) are therefore an interesting possibility but do not seem to be necessary for an explanation. Of course, if such large-scale asymmetries of the progenitor structure occurred, they could have important consequences for the SN explosion and the NS kick mechanism. Unquestionably, the availability of more realistic, three-dimensional progenitor models for core-collapse simulations would be highly desirable.


1

We adopt this term following a suggestion by Jeremiah Murphy.

2

Unfortunately, 3D progenitor models with self-consistent information about the asymmetries present in the stellar core before collapse are not available yet.

3

We repeat, however, that we do not use a light-bulb approximation, but in our simulations the outgoing neutrino luminosities are considerably modified by accretion contributions to the values imposed at the inner grid boundary and listed in Table 1.

4

It is important to note that there is no momentum flux across the outer grid boundary.

5

The neutrino-driven wind that develops after the onset of the explosion is essentially spherically symmetric in our simulations. Possible angular variations caused by PNS convection have low amplitudes and short wavelengths and thus do not carry any significant amount of momentum. Also the rate of angular momentum loss in the dilute outflow is typically very low.

6

It should be noted that the small number of species in the employed network is likely to be responsible for a considerable overestimation of the 44Ti production in our simulations. A better determination of the 44Ti yields will require more sophisticated network calculations.

7

The displayed quantity has the advantage that the bulk of the produced nickel becomes visible. Both the nickel density, which essentially follows the radial gradient of the mass density, and the nickel mass fraction, which does not necessarily trace regions with high integral nickel mass, are less suitable for this purpose.

8

We note in passing that Podsiadlowski et al. (2004) as well as van den Heuvel (2010) erroneously attributed the acceleration of the NS to anisotropic neutrino emission instead of asymmetric mass ejection.

Acknowledgments

We thank Rodrigo Fernández, Thierry Foglizzo, Jerome Guilet, and Jeremiah Murphy for useful discussions, and the latter in particular for suggesting the term “gravitational tug-boat mechanism”. This work was supported by the Deutsche Forschungsgemeinschaft through the Transregional Collaborative Research Center SFB/TR 7 “Gravitational Wave Astronomy” and the Cluster of Excellence EXC 153 “Origin and Structure of the Universe” (http://www.universe-cluster.de). Computer time at the Rechenzentrum Garching (RZG) is acknowledged.

References

All Tables

Table 1

Chosen values of parameters for all simulated models (see Sect. 2.5 for definitions).

Table 2

Explosion and NS properties for all models (see Sect. 3.1 for definitions).

Table 3

Hemispheric ejecta yields of helium and metals, and total nickel mass for the high-kick models W15-1 and W15-2, and the moderate-kick models L15-1 and L15-2.

All Figures

thumbnail Fig. 1

Axis-free Yin-Yang grid plotted on a spherical surface demonstrating symmetric geometry of the Yin (red) and Yang (blue) regions. Both grids are the low-latitude part of a spherical polar grid, and therefore contain grid cells which are almost equidistant in angular directions.

Open with DEXTER
In the text
thumbnail Fig. 2

Entropy-isosurfaces (left) of the SN shock (grey) and of the high-entropy bubbles (green), and entropy distribution in a cross-sectional plane (right) at t = 140,248,515 ms, and 1.3 s for model W15-2. The viewing direction is normal to the xz-plane. Length scales are given by the yardsticks below each image. Small Rayleigh-Taylor mushrooms start growing at  ~100 ms (top panels). The rising high-entropy bubbles merge and rearrange to form larger bubbles by the time the explosion sets in (second from top). The NS starts accelerating due to the asymmetric distribution of the ejecta, the acceleration reaching its maximum at  ~500 ms (third from top). At this epoch, the ejecta show a clear dipolar distribution with more dense, low-entropy material concentrated in the hemisphere containing the kick direction. The final NS kick direction has already found its orientation, and its projection onto the xz-plane is shown by the white arrows pointing to the lower right direction. The bottom panels show the entropy structure when the essentially spherically symmetric neutrino-driven wind has developed, visible as the green central region enclosed by the highly aspherical wind-termination shock (lower right panel).

Open with DEXTER
In the text
thumbnail Fig. 3

Entropy-isosurfaces (left) of the SN shock (grey) and of the high-entropy bubbles (green), ray-casting images of the density (middle), and entropy distribution in a cross-sectional plane (right) for models W15-1 (top) and L15-2 (bottom) at about 1.3 s and 1.4 s after core bounce. The viewing directions are normal to the plane spanned by the NS kick and spin vectors of each model, which also defines the plane for the entropy slices. The SN shock has an average radius of 13 000 km and 14 000 km for models W15-1 and L15-2, respectively (a length of 5000 km is given by yardsticks below the left images). The kick and spin directions are indicated by the white and black arrows, respectively, in the middle figures. The NS locations are marked by black crosses in the right plots. The images in the middle correspond roughly to the projections of the density distribution on the observational planes. Dilute bubble regions are light-colored in white and yellow, while dense clumps appear more intense in reddish and bluish hues. The purple circular areas around the NS represent the dense inner region of the essentially spherically symmetric neutrino-driven wind. The wind is visible in green in the right images and is bounded by the aspherical wind-termination shock. The NS is clearly displaced from the geometrical center of the expanding shock in the direction of the kick vector pointing to the lower left for model W15-1. The NS in model L15-2 has a much smaller kick velocity than in model W15-1, and does not show any clear displacement but remains located roughly at the center of the expanding ejecta shell.

Open with DEXTER
In the text
thumbnail Fig. 4

Long-time evolution of the NS velocity (top) and of the corresponding NS displacement (bottom) for models W15-6, W15-1, W15-2, L15-1, and L15-2. The displacement of the kicked NS from the coordinate center, although not taken into account in our simulations, can be estimated by time integration of the NS velocity. Within 3.3 s the NS in model W15-6 exceeds a velocity of 700 km s-1 and continues to be accelerated at a modest rate. Note also that the NS velocities in models L15-1 and L15-2 have nearly reached their terminal values after 3 s. In all cases the NS shift relative to the coordinate center is negligibly small before the explosion takes off, and even in model W15-6 with its most extreme NS kick it is always dwarfed by the diameter of the shock wave and postshock ejecta, which expand with typically  ~10 000 km s-1.

Open with DEXTER
In the text
thumbnail Fig. 5

Time evolution of the normalized spherical harmonic expansion coefficients for different (ℓ,m) modes computed for the shock surface Rs(θ,φ) (left) and for the surface mass density D(θ,φ) (Eq. (2); right) of models W15-2 (top) and L15-2 (bottom). The polar axis of the coordinate system chosen for the expansion is aligned with the final NS kick direction. Note the dominance of the (ℓ,m) = (1,0) mode for model W15-2 with the highest kick velocity, while this dominance is absent in the case of model L15-2.

Open with DEXTER
In the text
thumbnail Fig. 6

Left: time evolution of the NS acceleration in model W15-2. The acceleration computed by the assumption of momentum conservation, i.e., by the negative derivative of the total linear momentum of the gas surrounding the NS (red; denoted by “mom.consv.”), is compared with the value obtained by the integration of all hydrodynamical and gravitational forces that act on a sphere of radius 1.3Rns around the PNS (black; denoted as “total”). Also shown are the absolute values of the individual contributions of all terms in Eq. (8), namely the momentum transfers by downflows (light blue), outflows (green), anisotropic pressure forces (purple), and the gravitational attraction by the asymmetric ejecta (blue). The last contribution clearly dominate the evolution over long periods. Right: NS velocity components in all coordinate directions and total value of the NS velocity. The two lines in each panel, which correspond to our two methods of computing the NS recoil acceleration in the left panel, can hardly be distinguished from each other.

Open with DEXTER
In the text
thumbnail Fig. 7

Ejecta morphology of the high-kick models W15-2 (left) and W15-6 (right) at approximately 3.3 s after bounce. The NS (indicated by black spheres, enlarged by a factor of 20 for better visibility) are accelerated by the gravitational tug of the asymmetrically distributed, dense ejecta to velocities of 575 km s-1 and 704 km s-1, respectively (green arrows pointing to the right). The ejecta “clumps” with the highest densities, which are mainly responsible for the gravitational attraction, are visualized by the red-colored surfaces of constant density of 104 g cm-3. The small black dot left of the NS in both images marks the coordinate origin. The NS displacement relative to this location was computed by time integration of the NS velocity. After  ~3.3 s it is only 1260 km in W15-2 and 1440 km in W15-6. The neutrino-heated, asymmetrically distributed inner ejecta are bounded by the semi-transparent, beige surface from outside and the semi-transparent, light blue surface from inside. The former surface is defined by a constant 56Ni mass fraction of 5%, the latter surface conincides with the termination shock of the neutrino-driven wind. The outermost, semi-transparent green surface (cut open towards the observer for better visibility of the inner volume) corresponds to the SN shock wave. The average shock radius is approximately 36 000 km at the displayed time. A yardstick of 10 000 km along with a tripod indicating the coordinate directions are provided in the lower right corner of each plot.

Open with DEXTER
In the text
thumbnail Fig. 8

NS kick velocities and acceleration values (indicated by the length of the arrows attached to the symbols) versus explosion energies for all computed models after 1.1–1.4 s of simulated postbounce evolution. Models corresponding to different progenitors are shown by different symbols as defined in the figure.

Open with DEXTER
In the text
thumbnail Fig. 9

Neutrino flux asymmetry at  ~150, 300, 600 ms after bounce for model W15-3 as representative case. The 4π-maps show the relative variations ΔF/ ⟨ F ⟩  in the total (i.e., sum for ν and of all three neutrino flavors) neutrino energy flux in all directions. Hot spots (yellow and red) during the postbounce accretion phase (t ≲ 400 ms; top and middle panels) are connected to strong convective downflows in the neutrino-heated postshock layer. The typical scale of the corresponding high-emission areas is 20–30 degrees, and the amplitude of variations ranges around 10–30% during most of the time, though rare, single maxima can reach up to  ~70%. After the accretion has ended (t ~ 400–500 ms), the flux variation pattern is determined by the cell structure of the convection inside the proto-neutron star (i.e., below the neutrinosphere) and the variation amplitute decreases to a few percent (bottom panel). The number in the lower left corner of the plots (dℰ/dt) gives the total neutrino energy loss rate as flux integral over the whole emitting surface.

Open with DEXTER
In the text
thumbnail Fig. 10

Left: entropy distribution in a cross-sectional plane of model L15-1 at 1000 ms after bounce. The plane is chosen such that it contains the final vectors of the NS kick velocity induced by aspherical mass ejection (vns, long white arrow pointing to the right) and of the kick velocity caused by anisotropic neutrino emission (vns, short white arrow pointing towards the 7 o’clock position). The relative lengths of both arrows correspond to the absolute values of the velocities (vns ≈ 160 km s-1 and vns ≈ 5 km s-1), but the short arrow for vns is stretched by a factor of 5 for better visibility. A yardstick in the lower left corner of the figure indicates a spatial length scale of 2000 km. The contribution of anisotropic neutrino emission to the total NS kick velocity is clearly a small effect. Right: time evolution of the neutrino-induced kick velocity vns and of the angle αkν between vns and vns. The final velocities point roughly in opposite directions.

Open with DEXTER
In the text
thumbnail Fig. 11

Time evolution of the NS kick velocity (upper left panel) compared to the absolute value of the NS angular momentum (lower left panel), the polar and azimuthal angles (θspin, φspin) of the NS spin direction (right top and middle panels), and the relative angle αsk between spin and kick directions (right bottom panel) for five selected models. Two of these models represent high-spin cases (N20-2, N20-3), two are low-spin cases (W15-6, W15-7), and one is an intermediate case (L15-3) from our model set. The build-up of the angular momentum occurs in mainly a single impulsive episode during the convectively perturbed accretion phase of the NS. Together with its angles the angular momentum reaches its terminal value long before the NS acceleration has ended. The clustering of all five plotted models around a final θspin ≈ 90° is a mere coincidence, and the whole set of computed models exhibits the expected random distribution in polar-angle space.

Open with DEXTER
In the text
thumbnail Fig. 12

56Ni distribution for model W15-2 at 365 ms and 515 ms after bounce as produced by explosive burning in shock-heated ejecta. The integrated nickel mass per solid angle is plotted in the cross-sectional plane that contains the NS kick (indicated by the black arrows pointing to the right) and spin vectors. Red colored regions show that more nickel is formed in the direction opposite to the NS kick direction. The size of the displayed volume can be estimated from the yardsticks given in the lower left corner of each image.

Open with DEXTER
In the text
thumbnail Fig. 13

Shock deformation for the high-kick models W15-1 and W15-2 and for the moderate-kick models L15-1 and L15-2 at 1.3 s or 1.4 s after bounce, as given in the upper right corner of each map. These times are close to the evolution phases when the nickel production in the shock-heated ejecta takes place. The color scale is linear with black indicating minima and red maxima. The differences between maxima and minima relative to the average shock radii are 30%, 35%, 34%, and 29%, respectively. Red crosses indicate the NS kick directions while white plus signs mark the directions opposite to the kicks. The white lines separate the corresponding two hemispheres. The northern hemisphere is defined as centered around the NS kick vector.

Open with DEXTER
In the text
thumbnail Fig. 14

Integrated nickel masses per solid angle in all directions displayed in 4π maps for models W15-1, W15-2, L15-1, and L15-2 at postbounce times well after all nucleosynthesis processes have seized in our simulations. Red crosses indicate the NS kick direction while white plus signs mark the direction opposite to the kick. The corresponding hemispheres are separated by a white equator line. Note the differences between the high-kick W15 models and the moderate-kick L15 models. The former show a clear enhancement of the nickel production in the hemisphere opposite to the kick, whereas the latter eject nickel significantly more isotropically.

Open with DEXTER
In the text
thumbnail Fig. 15

Volumetric three-dimensional visualization of the nickel distribution for models W15-1, W15-2, L15-1, and L15-2 at the postbounce time given in the top left corner of each panel along with the model name. The semi-transparent isosurfaces correspond to a chosen value of the nickel mass per grid cell of 3 × 1026 g and are displayed at a stage well after all nucleosynthesis processes have seized in our simulations. The orange vectors represent the NS kick directions and are scaled by the corresponding NS velocities (96 km s-1 for the shortest arrow and 575 km s-1 for the longest one). The high-kick models W15-1 and W15-2 in the upper two panels exhibit a clear asymmetry with much more nickel being ejected in the hemisphere opposite to the kick direction. In contrast, the moderate-kick models L15-1 and L15-2 in the lower two panels exhibit a more isotropic distribution of the nickel, in particular no obvious hemispheric asymmetry between kick and anti-kick directions. While the radial distribution of the nickel may be strongly affected and changed by subsequent mixing instabilities that develop after the outgoing shock has passed the composition-shell interfaces of the progenitor star, the hemispheric differences in the nickel ejection will not be destroyed during the later supernova explosion.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.