A&A 457, 281-308 (2006)
DOI: 10.1051/0004-6361:20054654

Two-dimensional hydrodynamic core-collapse supernova simulations with spectral neutrino transport

II. Models for different progenitor stars

R. Buras^1,2 - H.-Th. Janka¹ - M. Rampp^1, - K. Kifonidis¹

1 - Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85741 Garching, Germany
2 - Max-Planck-Institut für Physik, Föhringer Ring 6, 80805 München, Germany

Received 7 December 2005 / Accepted 30 June 2006

Abstract
Spherically symmetric (1D) and two-dimensional (2D) supernova simulations for progenitor stars between 11$~M_\odot$ and 25$~M_\odot$ are presented, making use of the PROMETHEUS/VERTEX neutrino-hydrodynamics code, which employs a full spectral treatment of neutrino transport and neutrino-matter interactions with a variable Eddington factor closure of the ${\cal O}(v/c)$ moments equations of neutrino number, energy, and momentum. Multi-dimensional transport aspects are treated by the "ray-by-ray plus'' approximation described in Paper I. We discuss in detail the variation of the supernova evolution with the progenitor models, including one calculation for a 15$~M_\odot$ progenitor whose iron core is assumed to rotate rigidly with an angular frequency of 0.5 rad s^-1 before collapse. We also test the sensitivity of our 2D calculations to the angular grid resolution, the lateral wedge size of the computational domain, and to the perturbations which seed convective instabilities in the post-bounce core. In particular, we do not find any important differences depending on whether random perturbations are included already during core collapse or whether such perturbations are imposed on a 1D collapse model shortly after core bounce. Convection below the neutrinosphere sets in 30-40 ms after bounce at a density well above 10¹² g cm^-3in all 2D models, and encompasses a layer of growing mass as time goes on. It leads to a more extended proto-neutron star structure with reduced mean energies of the radiated neutrinos, but accelerated lepton number and energy loss and significantly higher muon and tau neutrino luminosities at times later than about 100 ms after bounce. While convection inside the nascent neutron star turns out to be insensitive to our variations of the angular cell and grid size, the convective activity in the neutrino-heated postshock layer gains more strength in better resolved models. We find that low (l = 1, 2) convective modes, which require the use of a full 180 degree grid and are excluded in simulations with smaller angular wedges, can qualitatively change the evolution of a model. In case of an $11.2~M_\odot$ star, the lowest-mass progenitor we investigate, a probably rather weak explosion by the convectively supported neutrino-heating mechanism develops after about 150 ms post-bounce evolution in a 2D simulation with 180 degrees, whereas the same model with 90 degree wedge fails to explode like all other models. This sensitivity demonstrates the proximity of our 2D calculations to the borderline between success and failure, and stresses the need to strive for simulations in 3D, ultimately without the constraints connected with the axis singularity of a polar coordinate grid.

Key words: supernovae: general - neutrinos - radiative transfer - hydrodynamics

1 Introduction

The mechanism by which massive stars explode is still unclear. State-of-the-art models with a spectral treatment of the neutrino transport by solving the Boltzmann equation or/and its moments equations agree in the finding that in spherical symmetry (1D) neither the prompt bounce-shock mechanism nor the delayed neutrino-driven mechanism lead to explosions for progenitors more massive than about 10$~M_\odot$(e.g., Rampp & Janka 2002; Liebendörfer et al. 2001,2004; Thompson et al. 2003; Sumiyoshi et al. 2005). Previous two-dimensional (2D) simulations (e.g., Herant et al. 1994; Burrows et al. 1995; Janka & Müller 1996; Fryer 1999; Fryer & Heger 2000) and three-dimensional (3D) models (Fryer & Warren 2004,2002) show the importance of convective overturn in the neutrino-heating layer behind the stalled supernova shock, which can enhance the energy transfer from neutrinos to the stellar matter and thus cause "convectively supported neutrino-driven explosions''. These multi-D models, however, employed radical simplifications of the treatment of neutrinos, mostly by grey diffusion or in a parametric way as heating terms. Concerns about the reliability of such approximations of crucial physics in studies of the supernova explosion mechanism were expressed by Mezzacappa et al. (1998).

Also the influence of convective activity inside the nascent neutron star, i.e. below the neutrinosphere, on the explosion mechanism has long been a matter of debate and requires further studies. The Livermore group (Wilson & Mayle 1993,1988) obtained explosions in their basically 1D models by assuming that so-called neutron-finger mixing instabilities exist in the newly formed neutron star, which accelerate the energy transport from the neutron star interior to the neutrinosphere. Thus the neutrino luminosities are boosted and the neutrino heating behind the supernova shock is enhanced. The analysis by Bruenn & Dineva (1996) and more recently by Bruenn et al. (2004), however, has demonstrated that neutrino diffusion leads to lepton number equilibration between perturbed fluid elements and their surroundings that is faster than assumed by Wilson & Mayle (1993,1988). Therefore neutron fingers are unlikely to occur in the supernova core. Bruenn et al. (2004) instead discovered a new mode of doubly-diffusive instability, which they termed "lepto-entropy fingers'' and which is also associated with neutrino-mediated thermal and lepton diffusion. The importance of this phenomenon during the early, critical phases of the explosion, however, was recently questioned by Dessart et al. (2005) because of its slow growth compared to Ledoux convection. The latter, in turn, was predicted to play a role in supernovae on grounds of 1D models of the neutrino cooling phase of nascent neutron stars. A Ledoux-type of convection was indeed found to be present during the first second after neutron star formation in 2D hydrodynamic simulations by Keil (1997), Keil et al. (1996), Janka & Keil (1998) and Janka et al. (2001). The latter simulations, however, considered only the proto-neutron star without self-consistently following its feedback with the environment of the supernova core. Moreover, a grey, flux-limited equilibrium "ray-by-ray'' diffusion code for the neutrino transport was used, with strong simplifications in the description of the opacities.

Only recently multi-dimensional simulations of stellar core collapse and post-bounce evolution with a spectral treatment of the neutrino transport have become possible (Buras et al. 2003,2006; Livne et al. 2004; Walder et al. 2005; Swesty & Myra 2005a,b; Burrows et al. 2006). Although these current approaches are the first steps of removing the severe deficiencies of the previous generation of multi-dimensional models, all of them still contain approximations of various, and different, aspects in the treatment of 2D transport. Swesty & Myra (2005a,b), for example, use a flux-limited diffusion description, an approximation also made by Walder et al. (2005) and Dessart et al. (2005), who in addition solve the transport for all neutrino energy groups independently. In contrast, Buras et al. (2006) have developed a "ray-by-ray plus'' approximation based on a variable Eddington factor solver for the coupled set of neutrino moments equations and Boltzmann equation, including a full coupling of the energy bins by neutrino reactions and by Doppler and gravitational redshift effects.

The approximations employed by the different groups are diverse and might hamper a detailed quantitative comparison of the results in the near future, and might constrain such efforts to a purely qualitative level. Eventually it will be necessary to test and possibly replace the current approximations by a more rigorous solution of the transport problem in the five- or six-dimensional phase space and in a relativistic framework, once the corresponding codes have become available and the necessary substantial increase of computer power has happened (Cardall 2004; Cardall et al. 2005).

Here we present results obtained with the multi-dimensional neutrino-hydrodynamics code MUDBATH, which is the "ray-by-ray plus'' implementation of the PROMETHEUS/VERTEX code described in detail in Buras et al. (2005, Paper I). In continuation of our previous work (Buras et al. 2003,2006), where also a broader introduction into the status of the field and its open questions is provided, we present here 1D simulations for nine different progenitor stars with masses between 11.2$~M_\odot$ and 25$~M_\odot$, and compare them with 2D simulations for three of these stars. The core collapse and post-bounce evolution of these models was followed until nearly 300 ms after shock formation.

Using a state-of-the-art treatment of spectral neutrino transport for hydrodynamical supernova simulations, the main goals of our work are:

• we compare 1D and 2D models in order to obtain quantitative information about the influence of convection below the neutrinosphere on the neutrino emission, the evolution, and the structure of the nascent neutron star. We analyse the influence of proto-neutron star convection on the conditions in the neutrino-heating layer behind the shock, and assess quantitatively the impact of convective activity in the postshock layer on the possibility for reviving the stalled shock and for getting a delayed explosion;
• we also study the differences of convection between a non-rotating and a rotating 15$~M_\odot$ model, whose iron core spins rigidly before collapse with a period of about 12 s, leading to an "extreme'' period of the settled neutron star of the order of 1 ms. In addition, we investigate the effects of low-mode (dipolar, l = 1, or quadrupolar, l = 2) hydrodynamic instabilities during the post-bounce evolution of an 11.2$~M_\odot$ star, comparing simulations with a full 180 $\hbox{$^\circ$ }$ polar grid and simulations which are contrained to a $\sim$90 $\hbox{$^\circ$ }$equatorial wedge (with periodic angular boundary conditions), thus preventing the development of such low modes in the pattern of the fluid flow;
• moreover, we perform tests for the influence of (i) the numerical resolution, in particular in the lateral direction of our 2D polar grid, (ii) of the chosen size of the angular wedge, and (iii) of the way in which we perturb our models to initiate the growth of convective instabilities, i.e., whether we follow a perturbed 2D model through core collapse and core bounce, or whether we map a 1D model to the 2D grid shortly after bounce, imposing random perturbations at that time.

The paper is structured in the following way. Main results of our 1D supernova simulations for the chosen set of progenitor models - whose basic properties are compared in Appendix A - will be discussed in Sect. 2, supplemented with more details in Appendix B. The 2D models will be presented in Sect. 3, with an analysis of the effects of convection in the forming neutron star in Sect. 3.1, a discussion of convection in the neutrino-heating layer in Sect. 3.2, a description of our full 180 $\hbox{$^\circ$ }$ model in Sect. 3.3, of the rotating model in Sect. 3.4, and of neutrino emission anisotropies in Sect. 3.5. Most of our 2D simulations without rotation were started from 1D collapse models only shortly after bounce, at which time small random perturbations were imposed to seed convective instabilities. Since the adequacy of such an approach may be disputed, we also performed simulations where the collapse phase was followed in two dimensions. This allowed us to investigate the growth of inhomogeneities during infall and to assess the possible influence of that on the growth of convection after bounce (Appendix E). A summary and conclusions will follow in Sect. 4. Appendix C contains a linear analysis of the structural changes of the proto-neutron star which can be expected as a consequence of convection below the neutrinosphere, and Appendix D introduces a simple mixing scheme by which we achieved to reproduce in 1D simulations the main effects of proto-neutron star convection as observed in our 2D models.

   
2 One-dimensional models

We have chosen a total of nine progenitors from different groups doing stellar evolution modeling. Details of the models can be found in Appendix A. These models are listed there in Table A.1 and cover a zero age main sequence mass (ZAMS mass) range from 11 to 25 $~\ensuremath{M_{\odot}}$. They represent various types of pre-collapse stellar structures.

The one-dimensional core-collapse simulations were performed with our 1D neutrino-hydrodynamics code VERTEX with spectral neutrino transport, using spherical coordinates and the physics described in detail in Buras et al. (2006). This includes a state-of-the-art treatment of neutrino interactions and an approximative description of general relativistic effects. The simulations described in the previous and in the present paper were done with the equation of state (EoS) of Lattimer & Swesty (1991), supplemented by a general lepton-photon-baryon EoS (ideal gases with Coulomb corrections) that extends to densities below those described by the Lattimer & Swesty EoS The numerical resolution used for hydrodynamics and neutrino transport was also specified in Buras et al. (2006).

The evolution can be separated into the phases of collapse, bounce, prompt shock propagation, neutrino burst, and accretion phase, which is in some cases accompanied by a transient shock expansion. None of our 1D simulations yields a prompt or delayed explosion.

The phases of core collapse, shock formation and propagation, and  $\nu _{{\rm e}}$ burst at shock breakout reveal only little differences between the progenitors because the core structure and properties are very similar or become very similar during collapse (for details, see Appendix B and also Liebendörfer et al. 2002). At the moment when the $\nu _{{\rm e}}$ burst is emitted the shock has reached an enclosed mass of $\sim$ $1.0~\ensuremath{M_{\odot}}$outside of which differences in the progenitor structure become larger. Therefore the mass infall rate through the shock begins to differ between the models. At the same time the postshock velocities become negative and the shock has converted to a stalled accretion shock. Since the rate of mass accretion by the shock is very high shortly after bounce, the shock is still pushed outward for some time due to matter accumulating behind it. This "passive'' expansion of the shock comes to an end when the mass accretion rate has dropped sufficiently and the neutrino cooling and settling of the accreted matter withdraw support from the shock. The conditions during this phase change only slowly. The situation can therefore roughly be characterized by steady-state conditions that depend on a number of parameters governing the accretion phase, i.e. the mass accretion rate through the shock, $\partial_t M_{\rm sh}$, the mass $M_{\rm PNS}$ and radius $r_{\rm PNS}$ of the proto-neutron star (PNS), and the neutrino luminosity $L_\nu$ (for a detailed discussion, see Appendix B). In recent analytic studies (Janka 2001; Arcones Segovia 2003) hydrostatic solutions are discussed which fairly well describe the nearly stationary accretion of the stalled shock as seen in the numerical simulations at later post-bounce times, provided that the model does not encounter quick changes such as a sudden drop of $\partial_t M_{\rm sh}$, or the onset of the explosion.

All 1D models fail to explode and the shock retreats to <100 km towards the end of our simulations (see top panel in Fig. B.4 in Appendix B). What is the effect of neutrino heating behind the shock? To discuss this question we consider two timescales, following Janka & Keil (1998), Janka et al. (2001), and Thompson et al. (2005). The advection timescale of matter falling inward from the shock to the gain radius is

$$\tau_{\rm adv}(t) = -\int_{r_{\rm gain}(t)}^{r_{\rm sh}(t)} \frac{1}{v_r(r,t)} {\rm d}r~,$$ (1)

where the gain radius $r_{\rm gain}$ is defined as the innermost radial position where neutrino heating dominates neutrino cooling, and v_r(r,t) is the radial velocity (which is negative when matter is accreted by the PNS). The heating timescale is

$$\tau_{\rm heat}(t) = \frac{4\pi\int_{r_{\rm gain}(t)}^{r_{\rm... ...\int_{r_{\rm gain}(t)}^{r_{\rm sh}(t)} Q(r,t) r^2 {\rm d}r~},$$ (2)

where Q is the net neutrino heating rate per unit of volume and

$$\varepsilon_{\rm bind}^{\rm shell}(r,t) \equiv\varepsilon(r,t) + \Phi_{\rm 1D}^{\rm enclosed}(r,t),$$ (3)

is the so-called local specific binding energy already introduced in Buras et al. (2006). Here, $\varepsilon=e_{\rm int} + \frac{1}{2}v^2$ is the specific energy, $e_{\rm int}$ is the internal energy per unit of mass, and v is the absolute value of the fluid velocity. The gravitational potential  $\Phi_{\rm 1D}^{\rm enclosed}(r,t)$ is calculated taking into account only the mass enclosed by the radius r.

While the advection timescale represents the time matter spends in the gain layer (between shock and gain radius), the heating timescale measures the time needed for neutrino heating to deposit an energy equivalent to the binding energy of the matter. Clearly, heating is of no importance as long as $\tau_{\rm adv} \ll \tau_{\rm heat}$. Thus, for obtaining a neutrino-driven explosion, the condition $\tau_{\rm heat} \la \tau_{\rm adv}$ must hold for longer than a time interval $\tau_{\rm heat}$ (see the discussion in Janka & Keil 1998; Janka et al. 2001; Janka 2001; and Thompson et al. 2005). Note that this condition is not necessarily sufficient for an explosion but it tells one when a visible shock expansion can be expected.

Looking at the timescales and their ratio (Figs. 1, 2), one recognizes that most models have long heating timescales and thus neutrino heating is inefficient in causing shock expansion. The ratio $\tau _{\rm adv} / \tau _{\rm heat}$ is always less than 1/2 except in the two low-mass models, s11.2 and n13, where neutrino heating is stronger during a short period of time.

	Figure 1: The advection timescales as defined in Eq. (1) versus post-bounce time. The lines are smoothed over time intervals of 5 ms.
Open with DEXTER

	Figure 2: The timescale ratios, $\tau _{\rm adv} / \tau _{\rm heat}$, as functions of time. The lines are smoothed over time intervals of 5 ms.
Open with DEXTER

The ratio becomes largest during the transient shock expansion associated with the times when composition interfaces reach the shock and the entropy makes a jump. This phenomenon, which is most extreme at the edges of the small iron cores in Model n13 at  $t_{\rm pb}\sim30$ ms and in Model s11.2 around $t_{\rm pb}\sim100$ ms, but in a weaker form is also present in Models s1b, s20.0, and s15s7b2 at around $t_{\rm pb}=120$ ms, 135 ms, and 170 ms, respectively, occurs because the density and therefore the mass accretion rate of the infalling matter drops at the interfaces of shells of different composition (Appendix A). The sudden decrease of the ram pressure allows for a transient shock expansion until the shock finds a new equilibrium position at a larger radius.

With a larger shock radius the postshock velocities are lower and the advection timescale increases. More efficient neutrino heating further strengthens the shock, the shock can therefore expand to even larger radii. This behaviour is obtained in Model n13, and especially in Model s11.2, where $\tau _{\rm adv} / \tau _{\rm heat}$ is close to unity for a period of about 20 ms. This, however, is shorter than the heating timescale  $\tau_{\rm heat}\simeq 30$ ms. The pronounced growth of the shock radius and of the advection timescale for Model s11.2, which results in the large local maxima of the corresponding curves in Figs. 1, 2, and B.4, is produced by a continuous strong decrease of the mass accretion rate, see Fig. B.4. After 113 ms post bounce this phase is over and the mass accretion rate continues to decline less rapidly. The shock then retreats quickly and finds a new quasi-stationary radius, where, however, neutrino heating becomes less efficient again.

In Models s15s7b2, s20.0, and s1b the drop of $\partial_t M_{\rm sh}$ leads also to shock expansion, but the effect is not strong enough to change  $\tau _{{\rm adv}}$ significantly; also, in these models the composition interface reaches the shock so late that the shock is already deep in the gravitational potential well and the postshock velocities are considerably higher than in Models n13 and s11.2. Therefore  $\tau _{\rm adv} / \tau _{\rm heat}$ remains well below unity.

In summary, we find that all our 1D simulations evolve both qualitatively and quantitatively in a similar way. In spite of maximum shock radii around 130-150 km the models do not reveal explosions. Only in the two lightest progenitor models, s11.2 and n13, the drop in  $\partial_t M_{\rm sh}$ at the composition shell interfaces is sufficiently steep and large and happens sufficiently early to allow the shocks to reach a radius of 170 km. Nevertheless, the models are far from producing explosions because the advection timescales remain always shorter than the timescales for neutrino heating, and the phases where the ratio of both approaches unity are much shorter than the heating timescale itself. Therefore neutrino heating is not strong enough to drive an explosion.

Table 1: Parameters of computed 2D models for progenitor stars with different masses. $\Omega _{{\rm i}}$ is the angular velocity of the Fe-core, which is assumed to rotate uniformly, prior to collapse, $\theta _0$ and $\theta _1$ are the polar angles of the lateral grid boundaries, and $N_{\theta }$ is the number of grid points in the lateral direction. $t_{\rm 2D}$ denotes the time (relative to the bounce time) when the simulation was started or continued in 2D.

The one-dimensional models analyzed so far have one significant shortcoming: they do not take into account hydrodynamic instabilities in the stellar core. Convection, especially the so-called hot bubble (HB) convection in the gain layer, has been seen to strengthen the shock in previous multi-dimensional supernova simulations. We can analyze our 1D models for the existence of Ledoux-unstable regions. For this purpose we introduce the new variable (defined as parameter $\chi$ in Foglizzo et al. 2005),

$$n_{\rm grow} (t) \equiv \int_{r_{\rm gain}(t)}^{r_{\rm sh}(t... ...ga_{\rm BV}(r',t)\right]_{>0}~ \frac{{\rm d}r'~}{v_r(r',t)} ,$$ (4)

where v_r is the radial component of the velocity and the Brunt-Väisälä frequency is defined as

$$\omega_{\rm BV}(r,t) \equiv {\rm sgn}\left(C_{\rm L}(r',t)\r... ...,t)} ~ \frac{{\rm d}\phi~ (r',t)}{{\rm d}r~} \right\vert ~~},$$ (5)

with ${\rm d}\phi~/{\rm d}r~$ being the local gravitational acceleration. $C_{\rm L}$ is the Ledoux-criterion, which is given by

$$C_{\rm L} = \left(\frac{\partial \rho}{\partial s}\right)_{{Y... ...e}}}\right)_{s,p} \frac{{\rm d}{Y_{\rm e}}~}{{\rm d}r~} \cdot$$ (6)

It predicts instability in static layers if $C_{\rm L}>0$. The Brunt-Väisälä frequency denotes the growth rate of fluctuations, if it is positive (instability), and the negative of the oscillation frequency of stable modes, if it is negative. In the gain layer between shock and gain radius, however, the gas is falling inward and the instability condition for Ledoux convection is not just given by $C_{\rm L}>0$(Foglizzo et al. 2005). Here the parameter  $n_{\rm grow}$ is of crucial importance and represents the number of e-foldings which short-wavelength perturbations will transiently be amplified during their advection from the shock to the gain radius. Since advection has a stabilizing influence, the threshold for convective instability in the gain layer is found to be $n_{\rm grow}\sim3$, and the growth of modes of lowest order (i.e., longest wavelengths) becomes possible only when $n_{\rm grow}\ga 5{-}7$(Foglizzo et al. 2005). Figure 3 shows that only some of our models get close to the critical value of  $n_{\rm grow}$ for which Ledoux convection can be expected according to a linear perturbation analysis. In case of Models s11.2 and n13, the value of $n_{\rm grow}$ clearly exceeds the critical threshold because in these models high entropy jumps at shell interfaces cause the shock to reach particularly large maximum radii (see Fig. B.4). This means that the postshock velocities entering the denominator in the integrand of Eq. (4) become smaller than in other models. It should be noted, however, that the discussion in Foglizzo et al. (2005) applies exactly only when the initial perturbations are very small. In case of our relatively large initial perturbations of order $1\%$, the fluctuations can grow to the non-linear regime - which is not accessible to the stability analysis - already for smaller values of  $n_{\rm grow}$.

Figure 3: Number of e-foldings that the amplitude of perturbations is estimated to grow during advection of the flow from the shock to the gain radius in our 1D models (cf. Eq. (4)). The lines are smoothed over time intervals of 5 ms. For Model n13, $n_{\rm grow}$ starts with values around 20 at $t_{\rm pb}=25$ ms, and begins to decrease at $t_{\rm pb}=35$ ms. Note that the entropy profile was smoothed before calculating the Brunt-Väisälä frequency, see the discussion in Buras et al. (2006, Sect. 2.5). The dotted line marks the instability criterion in advective flows according to Foglizzo et al. (2005), and corresponds to an amplification factor of about 20.

Open with DEXTER

   
3 Two-dimensional models

For our two-dimensional studies with assumed azimuthal symmetry we used the numerical code, input physics, and resolution as described in Buras et al. (2006) and specified in Table 1. The innermost core of 1.7 km radius was computed in spherical symmetry. We have chosen three representative progenitors: s11.2, which shows favorable conditions for developing strong Ledoux convection in the neutrino-heated postshock layer, as well as the two less promising progenitors, s15s7b2 and s20.0. A total of seven simulations were performed in 2D, see Table 1. Most simulations were started from a 1D model around 7 ms after bounce, at which time the radial velocity was randomly perturbed with an amplitude of $\pm$1%. For each progenitor we calculated a model with low angular resolution ( $2.7\hbox{$^\circ$ }$) with a $86.4\hbox{$^\circ$ }$ lateral wedge around the equatorial plane of the polar grid. In addition, we calculated high-resolution ( $1.35\hbox{$^\circ$ }$) versions for the two lighter progenitors s11.2 and s15s7b2. The corresponding 15$~M_\odot$ 2D  simulation was started at the onset of core collapse with an initial random density perturbation (in Model s15_64_p) of $\pm$2% in order to address the question whether the onset of convection after bounce changes when nonradial perturbations in the stellar core are followed through infall instead of being imposed shortly after bounce on a 1D collapse model in mapping the latter onto a 2D grid. Another simulation for the $11.2~\ensuremath{M_{\odot}}$progenitor, Model s112_128_f, with a resolution of  $1.41\hbox{$^\circ$ }$, was performed with a full  $180\hbox{$^\circ$ }$ grid. Finally, one simulation with the 15  $\ensuremath{M_{\odot}}$ progenitor, Model s15_64_r, included rotation. The angular frequency at the onset of core collapse was assumed to be $\Omega_{\rm i}=0.5 {\rm ~rad~s}^{-1}$ and constant in the Fe and Si core, and decreasing (spherically symmetrically) like r^-3/2 outside of 1750 km (1.43  $\ensuremath{M_{\odot}}$). This choice of the rotation rate and rotational profile was motivated by results for pre-collapse stellar cores of stars whose evolution is followed including the angular momentum transport by magnetic fields (Heger et al. 2005). Our choice of  $\Omega _{{\rm i}}$ in the iron core is roughly ten times lower than the core rotation of non-magnetic stars (Woosley et al. 2002) and about ten times larger than Heger et al.'s rotation rates of magnetized stars. It intends to maximize the effects of rotation during core collapse under the constraints that (a) the initial star can well be considered as spherically symmetric, and that (b) the assumed rotation, which is superimposed on a 1D stellar model, does not imply significant deviations from the hydrostatic equilibrium and gravitationally bound state of the progenitor model. In order to fulfill both requirements we limit the rate of rotation such that the influence of centrifugal forces is very small prior to collapse: $F_{\rm cent}/F_{\rm grav}<1\%$ everywhere for the initial ratio of the centrifugal force to the gravitational force. The density of the stellar model was perturbed by $\pm$1% and the simulation was carried out with a lateral resolution of  $1.41\hbox{$^\circ$ }$ in a wedge from the polar axis to the equator (i.e. besides axial symmetry also equatorial symmetry was assumed).

In Sect. 3.1 we shall focus on convection below the neutrinosphere of the proto-neutron star ("PNS convection''). We will try to classify this convection and will describe its effects on the evolution of the supernova and on the neutrino emission. To this end we will develop a mixing algorithm as described in Appendix D and will present a 1D simulation performed with it. We shall finish by discussing resolution and perturbation effects and differences between the simulations with different progenitors. In Sect. 3.2 a description of convective overturn in the neutrino-heating layer behind the shock ("hot bubble (HB) convection'') will follow. Again, we will first discuss differences compared to 1D models. Here, the resolution and size of the angular wedge play a much more important role and we need to elaborate on both aspects. Also the sensitivity of the onset of HB convection to the size of the seed perturbations will be discussed, and our results for the evolution of pre-collapse perturbations during the infall phase in Model s15_64_p can be found in Appendix E. Since differences of this model compared to Model s15_32 were minor, the results of this simulation will be included in the plots following below mostly without any special discussion. In Sect. 3.3 we shall describe our calculation with the full $180\hbox{$^\circ$ }$ grid, which allows low (l=1,2) modes to develop with important consequences. Section 3.4 will contain our results for the rotating $15~\ensuremath{M_{\odot}}$ model. We shall discuss the influence of rotation on the PNS, its neutrino emission, and convection according to the Solberg-Høiland-criterion (generalized to include the effects of ${Y_{\rm e}}$-gradients on convection). The role of centrifugal forces for the HB convection and shock propagation will also be addressed.

   
3.1 Proto-neutron star convection

Many phenomena inside the PNS, e.g. the conditions for convection in its interior, can mostly be discussed without taking into account the convective activity in the HB layer (but not inversely!). For several of the discussed models (s15_32, s20_32, s15_64_p) this is true, because they develop only weak HB convection with little impact on the shock propagation and on the PNS. Therefore nearly steady-state conditions prevail around the PNS, the mass in the gain layer changes only very slowly, and the mass accretion rate onto the PNS is approximately the same as the mass accretion rate through the shock. Moreover, the accretion flow onto the PNS is nearly laminar and isotropic.

Stability criterion.

In order to analyze the hydrostatic neutron star for convective instability we consider the Brunt-Väisälä frequency, Eq. (5). However, since neutrinos are strongly coupled with the dense plasma and close to equilibrium with the matter in the PNS, we generalize the Ledoux-criterion to a "Quasi-Ledoux criterion'', in which the effects of neutrino transport are approximately accounted for (see Wilson & Mayle 1993; Buras et al. 2006),

 

$$C_{\rm QL} \equiv \left(\frac{\partial\rho} {\partial s_{\rm tot}... ...m d}r~} - \beta_{\rm diff} \frac{{\rm d}{Y_{\rm lep}}~}{{\rm d}r~} \right) > 0,$$ (7)

for instability, where $s_{\rm tot}\equiv s+s_\nu$ is the entropy including the neutrino entropy, ${Y_{\rm lep}}\equiv {Y_{\rm e}}+ Y_\nu$ is the electron lepton number, and $\beta_{\rm diff}\sim1$ is a parameter which expresses the efficiency of lepton number exchange between buoyant fluid elements and their surroundings via neutrino diffusion. For our evaluation we choose $\beta_{\rm diff}=1$, following Keil (1997). The brackets $\left<\right>$ denote angular averages, whereas the last term, $\beta_{\rm diff}{\rm d}{Y_{\rm lep}}/{\rm d}r$, involves local evaluation because it describes the change of lepton number along the paths of buoyant fluid elements by the neutrino exchange between such fluid elements and their surroundings. Note that due to the nearly perfectly central gravitational acceleration in case of non-rotating neutron stars, only the radial components of the gradients of  $s_{{\rm tot}}$ and  ${Y_{\rm lep}}$ occurring in the general vector form of the Ledoux criterion (cf. Eq. (9)) are relevant (buoyancy acts only perpendicular to equipotential surfaces). Moreover, we simplified the analysis by performing it on the laterally averaged stellar background. The determination of stable and unstable layers on the basis of laterally averaged stellar quantities leads to essentially the same conclusions as the lateral average of the Ledoux criterion after evaluation with the 2D data.

Recently, Bruenn et al. (2004) presented an elaborate discussion of hydrodynamic instabilities in the PNS including the effects of neutrino diffusion (an extension of a previous analysis by Bruenn & Dineva 1996). They argue that local perturbations in the lepton number will be reflected in the neutrino phase space and thus cause a net neutrino diffusion which tries to wash out the perturbations, an effect which can be accounted for by a "response function''. Since neutrinos also carry entropy, the neutrino diffusion that smoothes the lepton number perturbations will create entropy perturbations. This effect is characterized by a "cross response function''. Of course, entropy perturbations will in an analogous way induce an equilibrating net neutrino diffusion which at the same time transports lepton number between the fluid elements and their surroundings, corresponding to another "response function'' and a "cross response function''. Bruenn et al. (2004) found in a numerical analysis that perturbation-induced neutrino lepton number transport by diffusion is considerably more rapid than thermal transport, and that the transport of lepton number reacts faster to entropy perturbations than to lepton number perturbations. For such a situation convective instability should set in for most stellar conditions, even when the fluid is Ledoux stable. In particular, Bruenn et al. (2004) describe two kinds of instabilities in the presence of neutrino diffusion: one instability occurs when the entropy difference between a displaced fluid element and its surroundings in a background with an entropy gradient results in a lepton fraction difference, which provides a driving force such that the induced perturbation grows. The buoyant rise of a perturbed fluid element, together with neutrino diffusion, will thus further increase the difference in entropy between the fluid elements and their surroundings and will create lepton number fluctuations from entropy perturbations, which continue to drive buoyant motion. A second instability exists where the neutrino diffusion creates an "overstable'' situation, i.e. where the effect of neutrino diffusion will drive a perturbed fluid element back to its original position, but to such an extent that the fluid element overshoots and thus oscillates around its original position with increasing amplitude. Bruenn et al. (2004) call these two doubly-diffusive instabilities "lepto-entropy finger'' (LEF) convection and "lepto-entropy semiconvection'' (LESC), respectively. They also distinguish Ledoux convection. However, Ledoux and LEF convection are closely related (LEF convection is an extension of Ledoux convection).

Bruenn et al.'s analysis of stellar profiles shows that Ledoux/LEF convection should appear in an extended region of the PNS from around 15 km to the neutrinospheres, whereas LESC should be visible deeper in the core. We suppose that our stellar profiles might yield qualitatively similar results if one applied their analysis. Our two-dimensional simulations, however, do not confirm their predictions in detail. For example, Model s15_32 shows convective instability between 17 km and 30 km at $t_{\rm pb}=48$ ms, see Fig. 5a. We interpret this as Ledoux convection and determine growth rates of the Ledoux instability (cf. Figs. 6 and 9) that are typically more than a factor of 10 higher than those of LEF found by Bruenn et al. (2004). Our hydrodynamic simulations reveal a behavior in agreement with our larger rate estimates. Moreover, the region showing convective activity does not extend to the neutrinospheres, which are beyond 60 km at that time (Fig. 11), and also at later times the upper boundary of the convective layer stays well below the neutrinospheres (Fig. 5b).

We suspect that the discussion by Bruenn et al. (2004) is not fully applicable up to the neutrinospheres. It is based on the assumptions that the mean free path (mfp) of neutrinos is much smaller than the size of the displaced fluid elements and shorter than local gradients in the background medium so that neutrinos diffuse and stay close to local equilibrium. In contrast, Fig. 4 reveals that the neutrino mfp becomes large near the upper edge of the convective zone, in particular for muon and tau neutrinos. As a consequence, the neutrinos start streaming increasingly rapidly and the neutrino densities in this region (first for $\nu _\mu$ and $\nu_\tau$, then for  $\bar\nu_{{\rm e}}$ and finally for  $\nu _{{\rm e}}$) begin to deviate from local chemical equilibrium. The net rate of neutrino losses therefore grows strongly, leading to a quick rise of the neutrino luminosities. Despite their growing mfp, the production of neutrinos in the stellar medium is sufficiently strong to allow for efficient energy and lepton number drain. We therefore suspect that in the semitransparent layer below the neutrinosphere the accelerating radial transport causes enhanced neutrino release such that entropy and lepton number fluctuations are wiped out by neutrinos streaming off. This happens before local exchange between fluid elements and their surroundings can efficiently take place as assumed by Bruenn et al. (2004). LEF instabilities thus have no time to grow.

Figure 4: Top: transport mean free paths of muon and tau neutrinos and antineutrinos for different energies as functions of radius at 243 ms after bounce, about 200 ms after the onset of PNS convection in Model s15_32. Also shown is the local density scale height (bold solid line). The vertical line marks the radius of the outer edge of the convective layer in the PNS (cf. Fig. 5). Bottom: mean free paths of $\nu _{{\rm e}}$, $\bar\nu_{{\rm e}}$, and $\nu _x$ as functions of neutrino energy at a radius r = 25 km (marked by the vertical line in the upper panel). The arrows indicate the mean energies of the $\nu _{{\rm e}}$ flux (solid), $\bar\nu_{{\rm e}}$ flux (dashed), and $\nu _x$ flux (dotted, coinciding with the dashed arrow) at this position in the star.

Open with DEXTER

Figure 5: Snapshots of PNS convection in Model s15_32 at 48 ms a) and 243 ms b) after bounce. The upper left quadrants of each plot depict color-coded the absolute value of the matter velocity, the other three quadrants show for $\nu _{{\rm e}}$, $\bar\nu_{{\rm e}}$, and $\nu _x$ (clockwise) the ratio of the local neutrino flux to the angle-averaged flux as measured by an observer at infinity. The diagonal black lines mark the equatorial plane of the computational grid and the thick circular black lines denote the neutrinospheres (which have a radius larger than 60 km in figure  a)). The figures c) and  d) show at the same times the relative lateral variations of lepton number and entropy (including neutrino entropy), i.e.  $\delta X \equiv (X-\left <X\right >_{\vartheta })/\left <X\right >_{\vartheta }$ for a quantity X.

Open with DEXTER

LESC as predicted by Bruenn et al. (2004) to occur far inside the PNS is not visible in our models. However, these authors mention that the existence of LESC is very sensitive to the exact values of the response functions, which in turn depend on the details of the neutrino interactions. We suspect that the different description of neutrino-matter interactions in our simulations might prevent LESC. In any case, convection very far inside the PNS should have less influence on the shock dynamics than the PNS convection which we see in our models. Therefore we will ignore the possibility of LESC in the following discussion and consider its effects on the neutrino emission and the explosion mechanism as negligible.

Phenomenology.

In Figs. 6-8, we see that a convectively unstable layer begins to develop in the PNS about 30 ms after bounce. This happens in a region which initially is stable due to a positive entropy gradient, see Fig. 7. As neutrino diffusion carries entropy away more efficiently at larger radii where the optical depth is lower, the entropy profile flattens, and finally the entropy gradient turns negative and the PNS becomes Ledoux unstable. Ten ms later the perturbations have grown (with a growth rate of about $\omega_{\rm BV}\sim 0.5{-}1.0{\rm ~ms}^{-1}$, see Figs. 6 and 9) to become non-linear and strong PNS convection has set in. This can be seen in the large lateral velocities in Fig. 8 (dark shaded). Similar to Keil (1997), we find that the structure of the convective cells is initially that of rolls with angular sizes between 20 $\hbox{$^\circ$ }$ and 30 $\hbox{$^\circ$ }$ and radial extension between 10 and 15 km (Fig. 5a). The rolls are stable for about 5 ms, then decay and form again at different locations in the same layer. 200 ms later, the contraction of the PNS has reduced the radial size of the convective cells to 10 km. The overturn velocities are around $3\times 10^8~~\ensuremath{{\rm cm}~{\rm s}^{-1}}$ with peaks of up to $5\times 10^8~~\ensuremath{{\rm cm}~{\rm s}^{-1}}$. At later times the velocities decrease to values around $2\times 10^8~~\ensuremath{{\rm cm}~{\rm s}^{-1}}$.

	Figure 6: Brunt-Väisälä frequency (Eq. (5)), using Eq. (7) as stability criterion, evaluated for different 1D models at 20, 30, and 40 ms after bounce.
Open with DEXTER

	Figure 7: Lepton number and total entropy versus enclosed mass in the PNS for the 1D model s15s7b2 for different post-bounce times before the onset of PNS convection in the corresponding 2D simulation.
Open with DEXTER

Figure 8: Convective region in Model s15_32. The dark-shaded regions have lateral velocities above $7\times 10^7~~\ensuremath{{\rm cm}~{\rm s}^{-1}}$, the light-shaded regions fulfill the Quasi-Ledoux criterion, Eq. (7). The solid lines mark the radii up to 30 km in 5 km steps, and 50 km; the dashed lines indicate density contours for 1014, 1013, 1012, and $10^{11}~\ensuremath{{\rm g}~{\rm cm}^{-3}}$.

Open with DEXTER

Figure 9: Brunt-Väisälä frequency in the 1D simulation of the stellar progenitor s15s7b2 (thin) and the corresponding 2D simulation, Model s15_32 (thick), at 30 ms (dashed), 62 ms (solid), and 200 ms (dash-dotted) post-bounce. At 30 ms, the thin and thick dashed lines coincide. All lines are truncated at 50 km. Positive values indicate instability. For the 2D model the evaluation was performed with laterally averaged quantities.

Open with DEXTER

The region with large lateral velocities (>700 km s^-1, dark shaded in Fig. 8) is wider than the layer with instability according to Eq. (7). This has two reasons: First, rising or sinking fluid elements can over- and undershoot into the adjacent stable layers. Second, the periodic boundary conditions applied in some of our simulations allow rings of uniform, lateral velocity to occur. This artificial phenomenon is associated with matter that settles from the neutrino-heated convective postshock layer onto the PNS and has obtained large lateral velocity components by participating in the overturn motion in the gain layer. The rings show up in Fig. 8 as a layer with high lateral velocities that moves from the "hot-bubble'' (HB) convective zone to the convection zone in the PNS between 100 and 150 ms post-bounce. Looking at the velocity distribution in 2D snapshots of the simulation, Fig. 5b, the rings can be identified around r=30 km.

Throughout our simulations, PNS convection occurs exterior to an enclosed mass of $0.6~\ensuremath{M_{\odot}}$, see Fig. 8. This inner boundary changes only little during the simulations, whereas the outer boundary of the convective layer moves outward in mass as time goes on, following the ongoing accretion of matter on the PNS. Keil et al. (1996) found in their models that PNS convection develops in an initially narrow layer, but the inner edge of this layer moves continuously deeper into the PNS. Their models, however, were evolved until 1.3 s after bounce, and the inward motion of the lower boundary of the convective zone may take place at times not covered by our present calculations. Moreover, the velocity at which the convective layer digs deeper into the PNS must be expected to depend on the treatment of the neutrino transport, which was described very approximatively by grey, radial, flux-limited equilibrium diffusion by Keil et al. (1996) and Keil (1997), who also made quite radical approximations for the neutrino-matter interactions (see Appendix in Keil & Janka 1995).

We note here that a closer discussion of the first 200-300 ms of post-bounce evolution, comparing these older diffusion simulations with our present models, is not useful to convey our understanding particularly of the possible deficiencies of the simpler diffusion treatment. Such a comparison is hampered by other very important differences between the simulations. Firstly, the post-collapse models from which the calculations of Keil et al. (1996) and Keil (1997) were started, were provided by Bruenn and thus were computed with a neutrino transport treatment that was not compatible with the one used for the subsequent 2D evolution. This must be suspected to have caused transient effects during an early phase of the post-bounce evolution of undetermined length. Such transients are not present in our current, consistent models. Secondly, the simulations by Keil et al. (1996) and Keil (1997) considered a neutron star of 1.1-1.2$~M_\odot$ and thus much smaller in mass than the compact remnants of the models in this paper, which grow in mass particularly rapidly during the first hundreds of milliseconds after bounce when the mass accretion rates are still high (see Fig. B.4). And thirdly, neglecting the accretion layer outside of the neutron star and placing the outer grid boundary near the neutrinosphere by Keil et al. (1996) and Keil (1997) may also have caused differences of their results from the present, full simulations of the collapsing supernova core. We therefore refrain from attempting a closer comparison of these old and our new models beyond the level of qualitative statements about the existence, long-lasting presence, and basic structural features of a convective layer inside the nascent neutron star.

The convective flow transports energy and lepton number from deeper layers of the PNS closer to the neutrinospheres, thereby flattening the entropy gradient and the lepton number profile, see Fig. 10. However, at the end of the calculation the initially quite steep negative  ${Y_{\rm lep}}$ gradient has not disappeared completely. This suggests that the convective mass motions can not efficiently transport lepton number over large distances. The reason for this is the fact that buoyant mass elements exchange leptons easily with their local surroundings via neutrino diffusion. Therefore rising bubbles with initially large  ${Y_{\rm lep}}$ end up with low  ${Y_{\rm lep}}$ when they reach the outer stable layers of the PNS. For this to happen the timescale of ${Y_{\rm lep}}$-equilibration between buoyant bubbles and their surroundings via neutrino diffusion must be of the same order as the rise timescale of the bubbles. This situation corresponds to a value of  $\beta_{\rm diff}$ in Eq. (7) close to unity. If the rise time were much longer, the 1D and 2D profiles of ${Y_{\rm lep}}$ in the convective region would not differ, while if it were much shorter, the ${Y_{\rm lep}}$ profile would be flat like the entropy profile.

A second effect contributes in establishing the observed conditions in the convective PNS layer. In contrast to the energy loss from the stellar medium, which is small in the considered period of time compared to the huge heat capacity of the accretion layer, lepton number is radiated away very efficiently, at least initially. At early times, when  ${Y_{\rm lep}}$ is still large in the convective PNS layer, the number flux of  $\nu _{{\rm e}}$ is much larger than that of  $\bar\nu_{{\rm e}}$ and the lepton number loss proceeds very rapidly. At later times, however, it slows down because ${Y_{\rm lep}}$ has decreased and the release of $\nu _{{\rm e}}$ and $\bar\nu_{{\rm e}}$ has become more similar.

	Figure 10: Lepton number, total entropy, specific internal energy, and density profiles versus enclosed mass in the PNS for a sample of 1D (dotted) and 2D simulations of different progenitor stars, 200 ms after bounce (i.e., approximately 160 ms after the onset of PNS convection). For the 2D models angle-averaged quantities are plotted.
Open with DEXTER

	Figure 11: Radius of the electron neutrinosphere (as a measure of the PNS radius) for the 2D models and their corresponding 1D models. For Model s15_64_r we show the equatorial ( upper) and polar ( lower) neutrinospheric radii. For other 2D models, angle-averaged quantities are plotted.
Open with DEXTER

Effects on radial structure and neutrino emission of the PNS.

The redistribution of lepton number and energy gradually affects the structure of the PNS. The "drain'' region, where ${Y_{\rm lep}}$ and internal energy $e_{\rm int}$ are extracted, has densities close to the nuclear density ( $0.3{-}2\times 10^{14}~\ensuremath{{\rm g}~{\rm cm}^{-3}}$), whereas the "dump'' region, where lepton number and energy are deposited, has lower densities ( $10^{12}{-}10^{13}~\ensuremath{{\rm g}~{\rm cm}^{-3}}$), see Fig. 8. A detailed analysis presented in Appendix C reveals that the transport of lepton number from the drain region to the dump region will effectively lead to expansion. Also the transfer of energy from the drain region to the dump layer leads to an expansion of the PNS. Consequently the neutrinosphere radii in 2D simulations are larger than in the corresponding 1D models, see Fig. 11.

The changes in the radial profiles of energy and lepton number and the expansion of the PNS relative to 1D simulations have interesting consequences for the neutrino emission. Convection in the PNS alters the emission of neutrinos in and above the dump region and therefore the PNS loses energy and lepton number at a different rate than without convection. In this case the linear analysis (Appendix C) shows that an increased loss of lepton number in and above the dump region leads to a more extended PNS, while an increased loss of energy leads to a more compact PNS (Fig. C.1). As will be discussed in detail below, PNS convection causes enhanced lepton number release and initially reduced energy loss via neutrinos, both supporting the "expansion'' of the PNS relative to the 1D simulations. At times later than 150 ms post bounce, the energy loss in neutrinos is also enhanced relative to the 1D models, without however being able to compensate the inflative effects on the PNS structure caused by the convective redistribution of energy and lepton number and by the enhanced lepton losses form the neutrinospheric layer due to neutrino emission.

Immediately after the onset of PNS convection the increase of ${Y_{\rm lep}}$near the upper boundary of the convective region leads to a higher electron chemical potential $\mu_{\rm e}$ and therefore a higher electron neutrino chemical potential $\mu_{\nu_{\rm e}}= \mu_{\rm e} + \mu_{\rm p} - \mu_{\rm n}$. Thus the $\bar\nu_{{\rm e}}$ abundance (given by equilibrium conditions) decreases significantly in this layer. Although this happens in a region where the $\bar\nu_{{\rm e}}$luminosity has reached only 5-10% of its final value (see Fig. 5a) the decrease is sufficiently large to lower the total $\bar\nu_{{\rm e}}$ luminosity by several percent. When the PNS convection has developed to full power, $L_{\bar\nu_{\rm e}}$ can be lowered by up to 10% in comparison to the 1D models, see Fig. 12. The decrease of the $\bar\nu_{{\rm e}}$ equilibrium abundance also affects the heavy-lepton ("$\nu _x$'' for $\nu _\mu$, $\bar\nu_\mu$, $\nu_\tau$, or  $\bar\nu_\tau$) luminosity via the process ${\nu_{\rm e}}{\bar\nu_{\rm e}}\rightarrow {\nu_{x}}{\bar\nu_{x}}$, which is the dominant $\nu _x$ production process in the region where the reduction of the $\bar\nu_{{\rm e}}$ abundance happens. However, $L_{{\nu_{x}}}$ in that region has already attained 70% of its final value, and therefore the reduced $\nu _x$ production rate leads to a decrease of the radiated $\nu _x$ luminosity by at most 5%, typically less (Fig. 12).

Figure 12: Luminosities of electron neutrinos, $\nu _{{\rm e}}$, ( top) electron antineutrinos, $\bar\nu_{{\rm e}}$, ( middle), and one kind of heavy-lepton neutrinos, $\nu_\mu,~\bar\nu_\mu,~\nu_\tau,~\bar\nu_\tau$, (denoted by $\nu _{{\rm x}}$; bottom) for the 2D models and for their corresponding 1D counterparts, evaluated at a radius of 400 km for an observer at rest. Note that the luminosities of Model s15_32 were corrected for the differences arising from the use of a slightly different effective relativistic potential as described in the context of Table 1.

Open with DEXTER

After 100 ms post-bounce these moderate effects are overridden by the structural changes of the PNS, which lead to larger neutrinospheric radii than in 1D models, as well as the ongoing convective transport of energy into the region below the $\nu _x$-sphere. Moreover, the convective layer now extends to lower densities so that the region affected by PNS convection contributes now 30%, 30%, and 90% to the radiated  $\nu _{{\rm e}}$, $\bar\nu_{{\rm e}}$, and $\nu _x$ luminosities, respectively (see Fig. 5b), although this region has a transport optical depth (for a definition see Buras et al. 2006, Eq. (28)) in excess of 40, 14, and 10 for  $\nu _{{\rm e}}$, $\bar\nu_{{\rm e}}$, and $\nu _x$, respectively, and therefore is still well inside the neutrinosphere.

Larger neutrinospheric radii without the described convective inflow of energy would lead to lower luminosities as a consequence of an associated decrease of the neutrinospheric temperature $T_\nu$. This, for example, is seen in simulations with different EoSs where the PNS radius depends on the high-density EoS properties. A larger PNS radius correlates with lower luminosities (Janka et al. 2005). In contrast, in our 2D models the luminosities increase. We indeed find lower $T_\nu$, which result in lower mean neutrino energies $\left< \varepsilon_\nu \right>$ (defined by the ratio of energy to number flux), see Fig. 13. The difference can be up to 10% for all neutrino kinds after 200 ms of PNS convection. Because of the energy transport to the neutrinospheres by convection, however, this reduction in $T_\nu$ is much weaker than it would be in an adiabatically expanding layer. Apparently, the larger neutrinospheric radii and only slightly lower temperatures lead to a net increase of the luminosities relative to the 1D results, see Fig. 12. The effect is strongest for $\nu _x$, which decouple energetically from the medium already near the upper boundary of the convective layer (Fig. 5b); after 200 ms of convection, $L_{{\nu_{x}}}$ is almost 20% higher than in the 1D models. $L_{\nu_{\rm e}}$ increases only by a few percent, while $L_{\bar\nu_{\rm e}}$ is almost identical in 1D and 2D  models, which means that the higher electron chemical potentials and the effects associated with the convective energy transport and structural changes of the PNS nearly compensate each other.

Figure 13: Average energies of the radiated neutrinos (defined by the ratio of energy to number flux) for the 2D models and for the corresponding 1D models, evaluated at a radius of 400 km for an observer at rest. The lines are smoothed over time intervals of 5 ms. Note that the average neutrino energies of Model s15_32 were corrected for the differences arising from the slightly different effective relativistic gravitational potential as described in the context of Table 1.

Open with DEXTER

After the onset of PNS convection the 2D models deleptonize faster than their 1D counterparts (Fig. 14). The lepton number loss is enhanced after 250 ms of post-bounce evolution by typically 8-10% (compare Fig. 14 with Fig. B.5). The evolution of the total energy loss is more complex and is smaller than in the 1D simulations for the first $\sim$100-140 ms of reduced energy emission. Only afterwards the losses become stronger in the 2D cases. However, even 250 ms after bounce, the total energy loss is only 2-4% higher with PNS convection.

	Figure 14: Differences between the total lepton number ( top) and energy losses ( bottom) of the 2D models and their corresponding 1D models as functions of post-bounce time. Here, $\cal {L}$ is the total lepton number flux, and L is the total neutrino luminosity.
Open with DEXTER

Note that the neutrinospheric luminosities in 2D models relative to 1D models are increased by 25%, 15%, and 25% for  $\nu _{{\rm e}}$, $\bar\nu_{{\rm e}}$, and $\nu _x$ after 200 ms. The effect of the PNS convection therefore manifests itself more strongly at the neutrinospheres than farther outside, because the luminosity differences get reduced due to structural differences above the neutrinospheres and thus resulting differences in the neutrino absorption and emission outside of the neutrinospheres.

We have performed several 2D simulations with varied resolution (2.7 $\hbox{$^\circ$ }$ and 1.35 $\hbox{$^\circ$ }$ or 1.41 $\hbox{$^\circ$ }$) and different choices of the lateral grid (either a wedge with  $86.4\hbox{$^\circ$ }$ around the equator or a full 180 $\hbox{$^\circ$ }$ grid). The results show that these differences affect the behaviour of the PNS convection only little. The flattening of the ${Y_{\rm lep}}$ profile appears to be slightly more effective at higher resolution, see Fig. 10, but the difference is too small to be distinguishable in other observables. Thus, a resolution of 2.7 $\hbox{$^\circ$ }$ seems sufficient to simulate PNS convection, mainly because the convective cells have a size of 20 $\hbox{$^\circ$ }{-}30\hbox{$^\circ$ }$ and are therefore much larger than the size of the angular grid zones. Similarly, the convective cells are sufficiently small to fit several times into the angular wedge in each of our simulations. Boundary conditions are therefore found to have no visible influence, and our simulation with a full 180 $\hbox{$^\circ$ }$ grid revealed no important differences of PNS convection and of its consequences compared to the models with a 90 $\hbox{$^\circ$ }$ angular wedge.

The PNS convection also evolves in a similar way in our simulations with different progenitors. The relative differences compared to the 1D models are quite similar for all progenitors. Moderate differences, e.g. in the total energy loss, $\int L^{\rm 2D} {\rm d}t~ - \int L^{\rm 1D} {\rm d}t~$, between Models s11_32, s15_32, and s20_32, see Fig. 14, appear small considering the large differences in the mass accretion rates and post-bounce evolution of these stars. The evolution of the energy loss is clearly different and lower only in the rotating Model s15_64_r, where a significant amount of energy is stored in rotation instead of being radiated by neutrinos. Also different seed perturbations have little effect on the evolution of convection in the PNS. Large seeds as, e.g., in Model s15_64_p, lead to PNS convection which starts typically 10-20 ms earlier. This short time difference of PNS convection, however, has hardly any influence on the supernova evolution.

Effects on the gain layer.

The structural differences between 1D and 2D models, which are induced by PNS convection, also influence the layer between the neutrinospheres and the shock. We can discuss these influences on the basis of Models s15_32 and s20_32, in which HB convection remains rather weak so that the differences from 1D result mainly from effects associated with PNS convection. In order to understand how much PNS convection contributes to the differences between 1D and 2D models we have also performed a 1D simulation, Model s15_mix, for the progenitor s15s7b2, in which PNS convection was treated with the simple mixing scheme described in Appendix D. In Fig. 10 we see that this approximative treatment reproduces the transport of energy and lepton number by PNS convection rather well. Therefore Model s15_mix can be understood as a simulation with PNS convection but without HB convection. Comparing Models s15s7b2, s15_32, and s15_mix, we conclude that model properties which basically are independent of the evolution of the gain layer, such as neutrinosphere radii and neutrino luminosities, are sensitive to PNS convection but not to HB convection (Figs. 11, 12). The differences between the 1D and 2D models in the gain region, which are visible in the mass of the gain layer, $M_{\rm gl}$, the advection timescale, $\tau_{\rm adv}^\ast$, and the shock radius, $r_{\rm sh}$ (Figs. 17, 15), are initially mainly caused by the HB convection, despite the fact that it is weak. However, the changes due to PNS convection become gradually more and more important, and make the dominant influence at the end of the simulations. In the following we will only discuss the effects which result from PNS convection. Model s15_mix will therefore also be considered as a "2D model''.

Note that in two-dimensional models with downflows and buoyant rising bubbles in the convective, neutrino-heated layer, Eq. (1) is not applicable as definition of the advection timescale, and we therefore use the following definition of an "effective advection timescale'':

$$\tau_{\rm adv}^\ast (M_i) \equiv \tau_{\rm adv}^\ast (t_1) = t_2(M_i) - t_1(M_i) ,$$ (8)

where t₂ is defined by the condition $M(r=r_{\rm gain},t=t_2) = M_i$ and t₁ by the condition  $M(r=r_{\rm sh},t=t_1) = M_i$. This is the time difference between the moment when the shock encloses a mass M_i and the time when this same mass M_i is enclosed by the gain radius. The expression is evaluated for different values of M_i to monitor the time evolution of  $\tau_{\rm adv}^\ast$. In cases where this definition of the advection timescale is used, we also apply it to the corresponding 1D models for comparison. In the multi-dimensional situation $\tau_{\rm adv}^\ast$ is a measure of the average storage period of gas in the gain layer, providing a rough estimate for the time interval which accreted matter is exposed to neutrino heating, without however considering in detail the complex fluid flow in downdrafts and rising plumes.

	Figure 15: Maximum and minimum shock radii as functions of time for the 2D models of Table 1, compared to the shock radii of the corresponding 1D models (dotted lines).
Open with DEXTER

As explained above, PNS convection leads to larger radii of the neutrinospheres. The mass enclosed by the neutrinospheres, however, is smaller because of a less compact PNS. The structure of the PNS is changed such that this also applies for the gain radius and its enclosed mass, as can be seen for Models s15s7b2, s15_mix, and s15_32 at late times in Fig. 16. However, the mass which is enclosed by the standing shock is identical in 1D and 2D models because it depends on the total mass accreted by the shock and therefore on the infall region and progenitor structure, which do not differ between the 1D and 2D models (the shock radius has almost no influence on the enclosed mass because of the relatively low densities behind the shock). For these reasons the mass in the gain layer must be expected to be larger in 2D. The top panel of Fig. 17 confirms this conclusion. We find that the mass in the gain layer, $M_{\rm gl}$, can be more than a factor of two larger than in 1D. Inevitably, with larger gain radius and larger mass in the gain layer, the shock radius is also larger than in 1D, see Fig. 15.

	Figure 16: Minimum and maximum gain radii as functions of post-bounce time for our $15~\ensuremath{M_{\odot}}$ 2D simulations compared to the gain radius of the corresponding 1D models.
Open with DEXTER

Figure 17: First ( top) panel: mass in the gain layer. Second panel: advection timescale as defined by Eq. (8). Third panel: total heating rate in the gain layer. Fourth panel: ratio of advection timescale to heating timescale. All lines are smoothed over time intervals of 5 ms. Note that the evaluation of  $\tau _{{\rm adv}}$ fails after the onset of the expansion of the postshock layer in the exploding Model s112_128_f.

Open with DEXTER

Since the density at a given radius between the neutrinospheres and the shock is higher in a model with PNS  convection (this is a structural consequence of the larger $\nu _{{\rm e}}$-sphere radius, which is located at $\rho\simeq 10^{11}~\ensuremath{{\rm g}~{\rm cm}^{-3}}$), and the mass accretion rate  $\dot M (r)$ is approximately constant for different radii above the gain radius (a feature of nearly stationary conditions), and $\dot M$ is equal in 1D and 2D (a consequence of the conservation of the mass flow through the shock), the postshock velocities are (on average) smaller in 2D. This implies a larger effective advection timescale, see Fig. 17, second panel. The values of  $\tau_{\rm adv}^\ast$ turned out to be larger by up to a factor of three in the 2D simulations.

Astonishingly, we find that the total heating rate of the gain layer, $\delta_t E_{\rm gl}$, is almost identical in 1D and 2D (Fig. 17, third panel). Model s15_mix, however, reveals that this is a coincidence. A comparison of this model with Model s15s7b2 shows that PNS convection actually reduces $\delta_t E_{\rm gl}$ slightly. We will not try to comprehend this behaviour analytically here, we only mention that the lower  $\nu _{{\rm e}}$ and $\bar\nu_{{\rm e}}$ luminosities and the lower neutrino energies (which the absorption rates depend on) lead to weaker heating. Also the fact that the gain radius is larger may contribute to this decrease. We believe that the larger mass in the gain layer is of minor importance for $\delta_t E_{\rm gl}$, because this additional mass is located at larger radii where the heating rate is small.

Since $M_{\rm gl}$ increases as a consequence of PNS convection, also $\tau_{\rm heat} \propto M_{\rm gl}/\delta_t E_{\rm gl}$increases. The ratio $\tau _{\rm adv} / \tau _{\rm heat}$, however, changes much less (Fig. 17, fourth panel). This can be understood by the fact that both  $M_{\rm gl}$ and  $\tau _{{\rm adv}}$ increase due to PNS convection, as explained above. These changes partly compensate each other in the evaluation of the timescale ratio.

In summary, PNS convection has two important consequences: First, the emitted neutrinos have lower mean energies than in 1D models (up to 10% lower after 200 ms of PNS convection). The neutrino luminosities initially decrease due to the onset of PNS convection ( $L_{\bar\nu_{\rm e}}$ by about 10% and  $L_{{\nu_{x}}}$ by about 5%) and increase at later times ($t\ga 150$ ms post-bounce), to reach several percent higher values for $\nu _{{\rm e}}$ and even $\sim$20% higher values for $\nu _x$ after 200 ms of PNS convection. Second, PNS convection affects the evolution of the gain layer and shock due to a less compact PNS and thus larger neutrinospheric radii with slightly lower temperatures and significantly higher luminosities near the neutrinospheres. This leads to a different structure between neutrinospheres and shock even without HB convection.

	Figure 18: Postshock convection in Model s15_32. Panel  a) shows snapshots of the entropy (in $k_{\rm B}$/by) at four post-bounce times. Panel  b) shows the radial velocity at the same times (the color bar is in units of 1000 km s-1 and its upper end corresponds to 2500 km s-1). The equatorial plane of the polar grid is marked by the black diagonal lines.
Open with DEXTER

   
3.2 Hot-bubble convection

We now turn to the discussion of Ledoux-type convection in the gain layer behind the shock, the so-called "hot bubble'' (HB) convection. Convective instability starts to develop in this layer in Models s15_32 and s20_32 at about 30 ms after bounce (Fig. 8), and convective anisotropies become visible, e.g. in the velocity and entropy distributions, shortly afterwards (see Fig. 18). The radial structure of the star, the shock radius (Fig. 15), and thus the evolution of the models after bounce are, however, hardly affected. In this sense the convective activity in the hot-bubble region is "weak''. The "growth number''  $n_{\rm grow}$ (Eq. (4), Fig. 3) for the growth of perturbations in the advection flow between shock and gain radius is correspondingly low. Obviously, neutrino heating is not powerful enough to drive strong buoyancy against the inward motion of the advected fluid. The small shock radius is associated with high negative postshock velocities, thus damping the growth of convection (for a general discussion, see Foglizzo et al. 2005). Only during the transient shock expansion, which sets in when the entropy jump at the Si-SiO interface reaches the shock, does the advection timescale become sufficiently long to give rise to a phase of slightly stronger overturn motions.

	Figure 19: Postshock convection in Model s112_64. Panel  a) shows snapshots of the entropy (in $k_{\rm B}$/by) at four post-bounce times. Panel  b) shows the radial velocity (in 1000 km s-1) at the same times, with maximum values of up to 35 000 km s-1 (bright yellow). The equatorial plane of the polar grid is marked by the diagonal lines.
Open with DEXTER

	Figure 20: Mass in the gain layer with the local specific binding energy as defined in Eq. (3) (but normalized per nucleon instead of per unit of mass) above certain values. The results are shown for Models s112_32 and s112_64 for times later than 100 ms after core bounce.
Open with DEXTER

More powerful HB convection occurs in Model s112_64. As we can see from Fig. 3, the perturbation growth factor increases to values of more than five in simulations with the s112 progenitor when the Si-SiO interface passes the shock 90 ms after bounce. The subsequent expansion of the shock is convectively supported and much stronger than in the 1D model (see Fig. 15). Although there is initially no obvious morphological difference in the convective pattern before and right after the shock expansion (compare Fig. 19, upper left panels with upper right panels), the conditions in case of a larger shock radius become more favorable for a significant strengthening of the neutrino-driven convection, in contrast to the situation in the more massive progenitors, where the composition interface has a much smaller impact. The long advection timescale of 30-40 ms (Fig. 17, 2nd panel) leads to a ratio of advection to heating timescale close to unity (Fig. 17, fourth panel). The mass in the gain layer with specific energy above certain values increases (Fig. 20), the pressure increases, and the shock is pushed to larger radii. The shock is highly deformed and the average shock position therefore fluctuates strongly due to the vigorous convective activity.

We find that in particular large-scale modes of the flow pattern gain strength when the gain layer becomes radially more extended. Besides Ledoux-instability of the neutrino-heated layer, the growth of these modes can be supported by the vortical-acoustic cycle in accretion flows (cf. Foglizzo 2001,2002; and Foglizzo et al. 2005), which has recently been discovered at work during the accretion phase of the stalled supernova shock (Blondin et al. 2003; Ohnishi et al. 2006; Scheck et al. 2004; Blondin & Mezzacappa 2006; Scheck et al., in preparation). The typical wavelength of the most unstable convective modes turns out to be linked to the radial width of the convective shell (for details, see Foglizzo et al. 2005). In Model s112_64 the number of rising high-entropy plumes decreases from five at 50 ms to only two at 107 ms and even only one at 115 ms after bounce (the equatorial bubble in Fig. 19 at 107 ms merges with the other bubble shortly afterwards). Such big structures persist for about 20 ms before they collapse and new rising bubbles form. The pattern is very nonstationary. Also at 140 ms and 170 ms large single bubbles appear, which dissolve again after roughly 20 ms. After the third generation of large bubbles has disappeared, the outer radius of the convective layer shrinks because of shock contraction, and smaller convective cells form again.

The sequence of generations of big, floating bubbles produces an oscillatory behavior of the shock radius (Fig. 15) as a consequence of a kind of feedback cycle between neutrino heating, bubble expansion, shock expansion and overshooting, reduced heating, shock contraction, bubble compression and collapse, increased accretion and heating, new bubble growth, and again shock expansion. This is similar to the oscillatory phases of shock expansion and recession, which we observed for one of the 1D models discussed in Buras et al. (2006, Sect. 3.1.4). In the present 2D situation, however, there is not the same coupling between shock behavior and neutrino luminosities found in that 1D model, where the accretion flow and thus the neutrino luminosities and neutrino heating were quenched during phases of rapid shock expansion. In the multi-dimensional case, downflows of accreted matter around a rising plume can still feed the cooling layer near the neutrinosphere and keep the neutrino luminosities at a high level (see Fig. 12). Here the expansion of the dominant rising plume lowers the heating rate of the bubble material, because its gas moves quickly away from the region of strongest neutrino energy deposition. This means that the expanding plume cuts off its own energy supply. The expansion is, however, still too weak to push the shock outward against the ram pressure of the infalling gas of the progenitor star. After transiently overshooting a possible equilibrium position, the shock therefore turns around, causing the compression of the high-entropy bubble and thus the onset of rapid cooling of the hot bubble gas by neutrino emission. The bubble collapses, and cool gas is channelled from the shock to the gain radius, increasing again the mass in the region of strongest heating. The energy input produces a new rising plume, which drives the shock outward again. The crucial question is whether the energy in the gain layer and the amount of matter with high specific energy increases or decreases during such a feedback cycle. In the first case an explosion may be produced as it happened in the 1D model of Buras et al. (2006) (where it was favored by a manipulation of the neutrino transport). In the present case, however, the feedback cycle ends after three periods (these three periods are visible in Fig. 15 as well as in Fig. 20). The decisive quantity is the ratio of advection to heating timescale, the time-average of which should be larger than unity to obtain net energy gain in the heating layer. For Model s112_64 this is not the case (Fig. 17, bottom panel). Although Model s112_64 gets much closer to success than any of the other simulations because of its strong convection and the large timescale ratio of almost unity, the simulation finally still fails to produce an explosion.

Turning now to resolution studies, we can compare Model s112_64 with Model s112_32. The latter has only half the lateral resolution but is otherwise identical with respect to input physics and radial grid. Clearly, the shock does not expand as far as in the high-resolution model (Fig. 15) and convection does not become as vigorous. Although the tendency of one convective cell to form with one rising plume and one downflow can also be observed in the less resolved model, the convective activity appears weaker and less neutrino-heated matter is involved in the overturn motions (Fig. 20). One reason for this is the fact that downflows become too narrow to be well resolved in the low-resolution model. Their lateral width and stability are therefore numerically overestimated, thus channelling more mass into the cooling layer below the gain radius. Numerical viscosity has a visible influence on the calculation.

Although the resolution in Model s112_64 is probably still not sufficient, the crucial quantity, the timescale ratio, increases only insignificantly when the lateral resolution is improved from 2.7 degrees to 1.35 degrees. We therefore have doubts that a further increase of the resolution can bring the timescale ratio above unity and thus establish favorable conditions for an explosion. Anyway, the resolution in Model s112_64 appears to be the minimal resolution needed for reasonably converged simulations. An angular cell width of one degree would be preferable according to systematic resolution tests (L. Scheck, private communication).

Figure 21: Postshock convection in Model s112_128_f. Panels  a)-c) show snapshots of the entropy for six post-bounce times. Panels  d)-f) display the radial velocity at the same times with maximum values of up to 47 000 km s-1 (bright yellow). In Panel  d) also the convective activity below the neutrinosphere (at radii $r \protect\la 30~$km) is visible. It is characterized by small cells, which are very similar to those observed in all other 2D models, too. The polar axis of the spherical coordinate grid is directed horizontally, the "north pole'' is on the right side.

Open with DEXTER

The size of the seed perturbations for starting convection turns out not to have an important influence on the HB convection. Comparing Models s15_32 and s15_64_p we find that a larger seed (as in the latter model) leads to HB convection which reaches the nonlinear regime about 10 ms earlier. The subsequent evolution, however, is very similar, and the faster onset of convective activity does not have any noticeable long-time effects. The reason may be the fact that the other relevant timescales, in particular $\tau_{\rm heat}$ and  $\tau _{{\rm adv}}$, are typically longer than the perturbation-dependent differences of the growth of HB convection.

   
3.3 A simulation with a grid of 180 degrees

We now discuss Model s112_128_f. The model is different from Model s112_64 in some important technical aspects. The full star was simulated with a lateral grid of 180 $\hbox{$^\circ$ }$ instead of the $\sim$ $90\hbox{$^\circ$ }$ wedge around the equatorial plane, which we used for Model s112_64. Along the polar axis of the spherical grid the gas flow is reflected, in contrast to our choice of periodic conditions at the lateral boundaries of the 90 $\hbox{$^\circ$ }$ wedge. Due to the assumed symmetry in case of the 90 $\hbox{$^\circ$ }$ grid, artificial constraints are imposed on the fluid flow. This has a selective influence on the type of flow pattern which can develop in the neutrino-heated hot-bubble layer. Global asymmetries with a dominant contribution of low modes - i.e., flow with l=1 (dipolar) or l=2 (quadrupolar) character in terms of orders of an expansion in Legendre polynomials for the cosine of the polar angle - can, for example, only be seen when the star as a whole is simulated. The possibility of long-wavelength modes in the postshock convective zone of the 11.2$~M_\odot$ simulation is in fact suggested by the linear perturbation analysis of Foglizzo et al. (2005), because  $n_{\rm grow}$ of Eq. (4) reaches values around 6-7 in our 1D run (Fig. 3). Low-mode convection and global asymmetries of the flow morphology are indeed found to grow and to cause significant quantitative differences in the evolution of Model s112_128_f compared to Model s112_64. These differences are large enough to change the outcome of the simulation for the 11.2$~M_\odot$ progenitor even qualitatively.

Noticeable differences between both models occur after about 60 ms post-bounce evolution (see Fig. 21). Although the overall morphology of the flow pattern in the convective postshock layer looks still similar in both simulations, bubbles at the poles of Model s112_128_f extend to larger radii. This is probably a consequence of the numerical setup of the simulation, which assumes axial symmetry and a polar grid with an axis that is impenetrable for the fluid flow. Flow which converges towards the axis is directed either inward or outward. Therefore its motion and behavior are constrained by the existence of the polar grid axis. Moreover, polar features have smaller volumes compared to structures near the equator, which are treated as tori around the symmetry axis. This geometrical difference of polar and equatorial structures is known to lead to differences in the growth rate of perturbations. Laser experiments (Kane et al. 2000) as well as 3D simulations (Scheck 2005) suggest that the smaller axial bubbles can grow faster than equatorial tori and behave more like the mushrooms of Rayleigh-Taylor instabilities in the truely 3D case.

In the further evolution of Model s112_128_f the shock and convective layer develop a large deformation with dominant dipole and quadrupole modes, which become more and more prominent. Huge bubbles inflate alternatingly in both hemispheres, while downflows are present near the equatorial plane. These downflows are very nonstationary and flutter back and forth between the hemispheres. Eventually, 180 ms after bounce, the gain layer is completely dominated by the two polar bubbles and one or sometimes two downflows around the equator, separated by a small transient bubble in between. The shock in this model reaches a radius of about 600 km at the end of our simulation (Figs. 15 and 22) and is expanding then with a speed that is typical of exploding models ($\sim$10 000 km s^-1). Unfortunately we had to terminate the simulation at 225 ms after bounce before it was possible to deduce the final parameters of the beginning explosion. The simulation had to be stopped because of a lack of computer time and the small timesteps enforced by large and rapid fluctuations of physical variables in the region where the equatorial downflow penetrates deep into the neutrinospheric layer.

	Figure 22: Shock radii at the poles and at the equator versus post-bounce time for Model s112_128_f.
Open with DEXTER

The polar plumes can be weakened transiently because of the ram pressure of the infalling material ahead of the shock or because of a descreasing supply of neutrino-heated matter when the downflows feed the high-entropy lobe in the opposite hemisphere. But the polar bubbles never collapse and contraction phases are reversed by new, powerful waves of high-entropy, high-velocity matter expanding away from the equatorial plane. The snapshot at 225 ms in Fig. 21 shows such a phase for the hemisphere on the left side (the polar grid axis is oriented horizontally): above 400 km the matter is still expanding, but below it has started retreating due to the ram pressure exerted by the matter falling through the shock. However, from below 300 km a new plume with very high velocities is already reviving the expansion of the gas between 300 km and 400 km. The whole situation is very dynamical.

The dipolar expansion is therefore driven and powered by the flow of neutrino-heated gas that is continuously replenished near the gain radius by the equatorial downflow of accreted matter. The consequences for the dynamical and energetic evolution of Model s112_128_f can be easily inferred from inspecting global quantities. In Fig. 17, bottom panel, we see that the timescale ratio $\tau_{{\rm adv}}^\ast/\tau_{{\rm heat}}$reaches values above unity in this model. A crucial difference of Model s112_128_f compared to its 90 $\hbox{$^\circ$ }$ counterparts is the fact that the advection timescale does not decline again after its maximum as in the other models but remains essentially constant for $t\ga 100~$ms post bounce. After a period of the order of $\sim$ $\tau _{{\rm adv}}$ the timescale ratio starts rising monotonically to climb to nearly 1.4 until 180 ms, associated with an expansion of the average shock radius (Fig. 22). This means that the dipolar instability has a stronger influence than Ledoux convection in the other models and leads to an enhancement of the efficiency of neutrino energy deposition. The total energy transfer by neutrinos becomes sufficiently large (in fact rises at $t\ga 170~$ms after bounce; Fig. 17) so that some part of the matter in the gain layer becomes nominally unbound and strong expansion sets in. Figure 24 shows that the mass of the gas with total local specific energy above some limits increases continuously, and the mass of the matter with positive energy follows this trend. Therefore we detect a steep growth of the "explosion energy'' (Fig. 24, upper panel) after 180 ms. This moment coincides with the accelerated expansion and marks the onset of the explosion. In Fig. 12 we see that the neutrino luminosities do not decrease compared to the luminosities of the non-exploding Model s112_64. This suggests that in spite of the launch of the explosion neutrino heating will go on at a significant level and will deliver more energy to the ejecta. Therefore Model s112_128_f does not exhibit the disadvantageous situation of explosions in spherical symmetry where the start of rapid shock expansion quenches the accretion of the forming neutron star and thus leads to a significant reduction of the neutrino luminosity.

The development of large-scale anisotropies with l=1,2 modes was also seen in other 2D simulations (Blondin et al. 2003; Ohnishi et al. 2006; Scheck et al. 2004; Burrows et al. 2006) and occurred in some preliminary 3D simulations (Mezzacappa & Blondin 2003; L. Scheck, private communication), too. The morphology of our 180 $\hbox{$^\circ$ }$ model looks similar to the models published by Blondin et al. (2003), who studied the hydrodynamics of nonradial instabilities of accretion shocks without taking into account the effects of neutrino transport and neutrino heating. However, we do not find that the growth of "turbulent energy'' in the expanding layer behind the shock plays a significant role for the explosion of Model s112_128_f. In contrast to the simulations by Blondin et al. (2003), Model s112_128_f shows a saturation of the kinetic energy of the lateral gas motions between neutron star (or gain radius) and shock at times later than $\sim$140 ms after bounce on roughly the same level as in Model s112_64. We neither find specifically large values of $E_{\rm kin,lat}$ in case of the 180 $\hbox{$^\circ$ }$ simulation, nor do we see a continuous growth or a distinct increase of that quantity correlated with the onset of the explosion and the development of positive values for the explosion energy after 180 ms (Fig. 23). We therefore conclude that the explosion of Model s112_128_f is driven by neutrino energy deposition and not by the amplification of turbulent kinetic energy in the postshock layer.

Figure 23: Kinetic energy, as a function of post-bounce time, associated with the lateral velocities of the matter between gain radius and shock (solid) and between the $\nu _{{\rm e}}$ sphere and shock (dashed) for Models s112_64 and s112_128_f. The rapid increase in case of Model s112_64 at $t_{\rm pb}\protect\ga 180$ ms is a numerical effect associated with lateral fluid flows which are enabled by our use of periodic boundary conditions at the angular grid boundaries.

Open with DEXTER

Figure 24: The upper panel displays the "explosion energy'' of Model s112_128_f, defined by the volume integral of the "local specific binding energy'' $\varepsilon _{\rm bind}^{\rm shell}$ as given in Eq. (3), integrated over all zones where this energy is positive, i.e.  $E_{\rm expl} = \int(\varepsilon_{\rm bind}^{\rm shell} > 0) \rho {\rm d}V~$. The lower panel gives the mass in the gain layer with the local specific binding energy (per nucleon) above certain values. The evolution of these quantities as functions of post-bounce time is shown. Note the difference from Fig. 20.

Open with DEXTER

Dipolar asymmetries seem to be generic phenomenon in supernova simulations in which the neutrino heating is not strong enough to revive the stalled shock within a short time, i.e. on a timescale of the order of the advection timescale in the gain layer. However, the question has to be asked whether the morphology of the flow is enforced in the 2D simulations by the choice of the azimuthal symmetry and the coordinate grid with its polar axis. The existence of a preferred axis direction predetermines and might even support the growth of a dipolar asymmetry in this direction. In contrast, 3D simulations without a preferred axis might reveal a slower development of such asymmetries. It therefore needs to be demonstrated that freedom of motion in all directions does not prevent the instabilities with dominant l=1, m=0 contributions. In reply to such objections one might say that in case of only a small amount of rotation the polar axis will obtain a physical meaning and will act as a centrifugal barrier for the flow (if the axial component of the specific angular momentum of fluid elements is conserved). The direction for the development of dipolar deformation might therefore naturally be selected by the axis of (slow) rotation. More arguments in favor of the physical nature of low-mode, global asymmetries in the supernova core come from analytic investigations and linear analyses. On the one hand these suggest that stalled accretion shocks are subject to nonradial instabilities with highest growth rates for l=1 deformation (Foglizzo 2002; Thompson 2000). On the other hand it was found that the l=1 mode is the fastest growing mode in case of volume-filling thermal convection in a sphere (Chandrasekhar 1961). The latter situation might apply to the supernova core, provided the shock radius is much larger than the inner boundary of the neutrino-heated, convective layer. But the true supernova conditions differ in important aspects, e.g., boundary conditions and accretion flow in the postshock layer, which prevent a direct quantitative application of Chandrasekhar's analysis (Foglizzo et al. 2005). In any case, 2D simulations as well as analytic and linear analysis can only be taken as suggestive. The relevance of low-mode asymmetries during the phase of shock revival will ultimately have to be demonstrated by 3D models without the handicap of a coordinate singularity along the polar axis of a spherical or cylindrical grid.

Therefore, we have to cautiously interpret the results of our present simulation. The choice of the lateral grid size and of the boundary conditions has turned out to decide whether the dipolar instability is suppressed as in case of the 90 $\hbox{$^\circ$ }$ wedge with periodic boundary conditions, or whether it can develop (and possibly is supported or even enforced) as in case of the 180 $\hbox{$^\circ$ }$ simulation with axial symmetry and reflecting boundary conditions along the grid axis. While Model 112_64 with the 90 $\hbox{$^\circ$ }$ wedge fails marginally, Model s112_128_f yields an explosion for the considered $11.2~\ensuremath{M_{\odot}}$ progenitor star. Because of the relatively low mass in the gain layer, the explosion might remain rather weak, but our simulation had to be terminated too early to be finally conclusive in this point. The quantitative similarity in many aspects but qualitative difference in the outcome of the two simulations demonstrates how close Model 112_64 was already to an explosion. Conversely, it can also mean that small effects which weaken the growth of the dipolar instability might delay the onset of the explosion (and thus change quantities such as the initial mass cut and the explosion energy), or might even lead to a failure.

Figure 25: Radial profiles of the specific angular momentum jz, angular velocity $\Omega$, and ratio of centrifugal to gravitational force, $F_{\rm cent}/F_{\rm grav}$, in the equatorial plane versus  a) radius and  b) mass for Model s15_64_r at different times. Note that the "enclosed mass'' in multi-dimensional simulations is defined here as the mass within a spherical volume of specified radius. Quantities plotted versus enclosed mass therefore do not represent Lagrangian information. Time is normalized to bounce. Note that in the upper right panel the lines for times -175 ms, 0, and 10 ms lie on top of each other because of specific angular momentum conservation and negligibly small rotational deformation of the collapsing stellar core until later after bounce.

Open with DEXTER

The bottomline is that our 11.2$~M_\odot$ simulation lingers at the borderline between failure and explosion. But the effect which has triggered the success in case of Model s112_128_f, i.e. the large-scale non-radial modes of the fluid flow between shock and gain radius, can be treated only approximately in two-dimensional simulations. We can not exclude that the success is a result of an overestimation of this phenomenon and of its consequences. Three-dimensional (axis-free) simulations with reliable neutrino transport are needed to convincingly demonstrate that the neutrino-driven mechanism, supported by low-mode convection and accretion shock instability, is viable to explain supernova explosions of massive stars at least for some range of progenitor masses.

Figure 26: Snapshots of Model s15_64_r. The rotation axis is oriented vertically. Panels a), b) show the distributions of the specific angular momentum jz, of the angular frequency $\Omega$, and of the ratio of centrifugal to gravitational force for the moment of bounce and for a post-bounce time of 271.1 ms, respectively. Note the different scales which represent enclosed mass M (which corresponds to values of mass enclosed by spheres of chosen radii) and radius r. The black lines mark the shock. Panels c), d) depict quantities which are of interest for discussing the physics in the gain layer at the time of maximum shock expansion and at the end of the simulation, respectively. The thick black line marks the shock, the thin black lines indicate the gain radius and $\nu _{{\rm e}}$ sphere. Panels e), f) provide analogous information inside the PNS at the onset of convection and at the end of the simulation, respectively. Instead of the gas entropy and ${Y_{\rm e}}$, the total entropy including neutrino contributions and the total lepton number are plotted. The three black lines in panel f indicate from outside inward the neutrinospheres of $\nu _{{\rm e}}$, $\bar\nu_{{\rm e}}$, and $\nu _\mu$.

Open with DEXTER

   
3.4 A simulation with rotation

We investigate the effects of rotation on the supernova evolution with Model s15_64_r, for which we assume that the iron core rotates with a period of about 12 s prior to collapse. This is rather slow compared to the conditions considered in many of the core-collapse simulations performed by other groups (e.g., Ott et al. 2004; Kotake et al. 2004). The reasons for our choice of the initial rotation law were explained in detail at the beginning of Sect. 3 and in Müller et al. (2004). The angular frequency $\Omega$ of the Fe-core at the beginning is 0.5 rad s^-1 (Fig. 25) and increases due to angular momentum conservation during collapse to maximum values around 600-700 rad s^-1 in the homologous core shortly after bounce, see Figs. 25, 26a. Rotation generates a centrifugal force  $F_{\rm cent}$ which at this time is at most $6\times 10^{-3}$ of the gravitational force $F_{\rm grav}$(Figs. 25, 26a). The ratio of rotational to gravitational energy, $\beta_{\rm rot} = E_{\rm rot}/\vert E_{\rm grav}\vert$, grows from an initial value of less than 10^-3 to roughly 0.4% at bounce (Fig. 27). Because of the subsequent contraction of the PNS and the inflow of material with higher specific angular momentum (j_z), the ratio $F_{\rm cent}/F_{\rm grav}$ rises up to 25% near the equator (Fig. 25) and can reach values around 3 in the polar regions at a post-bounce time of 271.1 ms (Fig. 26b). The angular frequency increases up to maximum equatorial values around 2000 rad s^-1, and even much larger values close to the polar axis (Fig. 26b).

As a consequence, the forming neutron star develops an increasing degree of rotational flattening. The oblateness of the neutrinospheres is most pronounced during the later stages of the simulated post-bounce evolution (Fig. 26, panels c, d, f and Fig. 29), whereas the deeper interior of the PNS shows a much smaller deformation, and the convective layer at $r \la 25~$km remains nearly spherical (Fig. 26, panels e, f). At the end of our simulation the polar value of the radius of the $\nu _{{\rm e}}$-sphere is about 30 km, which is roughly the same as in non-rotating models, whereas the $\nu _{{\rm e}}$-sphere is at 40 km near the equatorial plane (Figs. 11, 26). Also the rotation-generated differences of the gain radius in different directions grow with time, although the contrast visible in Fig. 16 shows the combined effects of rotation and local variations due to the convective anisotropies in the postshock region (see Fig. 26). The angular differences of the gain radius become roughly a factor of two until the end of our simulation for Model s15_64_r and are thus somewhat larger than in the other models with strong convective overturn in the neutrino-heating layer (in particular Models s112_64 and s112_128_f). The radial structure of the rotating PNS differs most strongly from that of the corresponding non-rotating Model s15_32 near the equatorial plane, where the centrifugal force, $F_{\rm cent}$, points in the opposite direction of the gravitational force, $F_{\rm grav}$. Differences in the radial structure between the rotating and non-rotating PNSs are less significant at the poles.

Figure 27: Ratio of rotational energy to the gravitational energy, $\beta _{\rm rot}$, as a function of time for Model s15_64_r. The quantity is shown for the whole computational volume (solid line) and for spherical volumes with an enclosed mass of the Fe-Si core (1.43  , dashed line) and of the Fe core (1.28  , dotted line), respectively. Note that in the latter two cases the time evolution of $\beta _{\rm rot}$ does not provide Lagrangian information, because matter is redistributed by rotational deformation and by convection.

Open with DEXTER

Because of the rotational flattening the PNS as a whole is less compact than in the non-rotating case, i.e., the volume enclosed by the neutrinospheres is larger. The implications of a bigger average neutrinospheric radius for the global properties of the neutrino emission were discussed in detail in Sect. 3.1. Due to the lower temperatures at the more extended neutrinosphere, the mean energies of the radiated neutrinos are reduced compared to the corresponding 1D model, associated with a decrease of the neutrino luminosities. In non-rotating models, this reduction holds only for the first $\sim$150 ms after bounce and is compensated later. In fact it is over-compensated by the convective energy transport to the neutrinosphere, which leads to higher luminosities (in particular for muon and tau neutrinos) in 2D simulations at $t\ga 150~$ms after bounce (in spite of continuously lower mean neutrino energies; Figs. 12 and 13). In the rotating Model s15_64_r this "cooling effect'' of the more extended neutrinosphere is significantly stronger than in the corresponding non-rotating Model s15_32, causing even lower energies of the radiated neutrinos (Fig. 13). In fact, the neutrinospheric temperatures are reduced so much by rotational expansion that this effect cannot be overridden by the convective transport of energy to the neutrinosphere, in particular also because PNS convection tends to be slightly weaker in Model s15_64_r than in Model s15_32 (compare Figs. 26 with 5 and see the discussion following later). The neutrino luminosities in the rotating model stay therefore clearly below those of its non-rotating counterpart during all of the simulated evolution (Fig. 12), and the integral of the radiated energy is smaller than in the 1D case (Fig. 14, lower panel), signalling that energy is stored in rotation instead of being converted to neutrino emission. In contrast, the rotating model does not show a particular behavior with respect to the lepton number loss (Fig. 14, upper panel), because the electron neutrino and antineutrino luminosities are both reduced by rotational effects in very similar ways (Fig. 12). The emission anisotropies caused by the rotational deformation of the neutrinosphere and by convective effects will be addressed in Sect. 3.5.

In contrast to the oblate deformation of the PNS, the shock reveals a prolate shape during most of the computed evolution, with a radius typically 50 km larger at the pole than at the equator (Figs. 15, 26, panels c, d). The polar bulge of the shock is maintained by a big convective bubble which exists rather stably between polar angles of about 10 $\hbox{$^\circ$ }$ and about 45 $\hbox{$^\circ$ }$. After an initial, transient phase in which the convective activity grows and smaller convective cells merge to larger, volume-filling structures, this large pole-near bubble comes to exist besides one or two additional plumes closer to the equator. Although these bubbles are again strongly time-dependent as in the other 2D models, and phases of bubble contraction are followed by bubble reinflation, the morphology is very stable and the bubbles grow again essentially at the same places. The pattern seems to be determined and supported by the presence of an angular momentum gradient in the neutrino-heated layer (Fig. 26, panel b) and by the action of centrifugal and coriolis forces on the fluid motion.

Towards the end of our simulation (from about 180 ms until 280 ms post bounce) quasi-periodic large-amplitude pole-to-equator oscillations with a cycle time of 15-25 ms set in (Fig. 15) where phases with a larger pole-near bubble alternate with time intervals in which the convective plume near the equator is stronger. The transition between both extrema is characterized by a merging of the two convective cells into one big plume at intermediate latitudes. The shape of the shock changes back and forth between a pronounced prolate deformation and a more oblate shape. During all these time-dependent variations the maximum and minimum shock radii stay around 200 km and 150 km, respectively (Fig. 15). The average shock radius at $t\ga 150~$ms p.b. is 50-100% larger than in the non-rotating 2D models and does not decay until the end of our simulation at nearly 300 ms after bounce. Due to centrifugal effects the presence of angular momentum in the infalling matter has a stabilizing influence on the postshock flow and on the shock. Thus rotation ensures a more extended gain layer and supports strong convection, in contrast to the non-rotating 15$~M_\odot$ Models s15_32 and s15_64_p, where the retraction of the shock after its maximum expansion strongly damps the convective activity in the gain layer (Fig. 15; and Buras et al. 2006). The convective pattern in the postshock region therefore depends on the amount of angular momentum carried by the accreted matter. Since the gas falling onto the shock comes from larger and larger initial radii at later times, the angular momentum of the accretion flow increases continuously (Fig. 25). This may explain some of the evolution which we observe for the morphology and dynamical behavior of the convective post-shock layer.

The conditions for convective instability are affected by rotation in a way which is expressed by the so-called Solberg-Høiland criteria (see e.g. Tassoul 1978; Keil 1997). For the considered situation in the supernova core, the first criterion is a generalization of our Quasi-Ledoux criterion, combined with the Solberg-term which accounts for the stabilizing effect of a positive angular momentum gradient. With the parameter $\beta_{\rm diff}$ for neutrino diffusion in the Quasi-Ledoux criterion (Eq. (7)), the corresponding mode frequencies are

 

$$\omega_{\rm Hoi,QL} = \pm \left\{ \left\vert - \frac{1}{x^3} \fra... ..._{p,s_{\rm tot}} \vec{\nabla} {Y_{\rm lep}}\right) \right\vert \right\}^{1/2} ,$$ (9)

where instability is given when $\omega_{\rm Hoi,QL} > 0$. The sign of this expression is identical to the sign of the absolute value term under the square root, $\vec{a}_{\rm cent}$ is the acceleration due to centrifugal forces, $s_{{\rm tot}}$ is the sum of gas and neutrino entropies, and $x=r \sin{\vartheta}$ is the distance from the rotation axis. Here, the Quasi-Ledoux criterion was generalized further to a fully multi-dimensional form by introducing the scalar product of the vector quantities. In the gain layer, the adequate Solberg-Høiland criterion is obtained from Eq. (9) by ignoring the neutrino entropy and lepton number, i.e. by replacing $s_{{\rm tot}}$ by s and  ${Y_{\rm lep}}$ by  ${Y_{\rm e}}$, and setting $\beta_{\rm diff}\equiv 0$ when neutrino diffusion effects are irrelevant. The second Solberg-Høiland criterion for instability is given by

$$C_{\rm Hoi}' \equiv {}- \frac{ (\vec{\nabla} \phi - \vec{a}_{\rm cent})_z}{\rho} \lef... ...t\vert _{p,s_{\rm tot}} \frac{{\rm d}{Y_{\rm lep}}~}{{\rm d}z~} \right] \right.$$

$${} \left. - \frac{{\rm d}j_z^2~}{{\rm d}z~} \left[ \left.\frac{\p... ..._{p,s_{\rm tot}} \frac{{\rm d}{Y_{\rm lep}}~}{{\rm d}x~} \right] \right\} > 0~,$$ (10)

where $z=r \cos{\vartheta}$ is the distance from the equatorial plane and the suffix z denotes the z-component of vector quantities. We will not further consider this criterion in the following because convective instability requires at least one of the two criteria to be fulfilled. An evaluation of our model, however, reveals that the second criterion indicates instability only when the first criterion is also fulfilled.

	Figure 28: ${Y_{\rm lep}}$ profile at 200 ms post-bounce time for a 1D model, a 2D model, and our 2D model with rotation. The plot shows angle-averaged quantities for the 2D cases.
Open with DEXTER

Figure 29: a) Mass $M_{\rm NS}$ and total angular momentum $J_{\rm NS}$ of the "neutron star'' (defined by the matter with density above $10^{11}~\ensuremath{{\rm g}~{\rm cm}^{-3}}$) in Model s15_64_r. b) Moment of inertia $I_{\rm NS}$ and average angular frequency $\overline {\Omega }_{\rm NS} = I_{\rm NS}/J_{\rm NS}$ of the neutron star in Model s15_64_r. c) Average radius $R_{\rm NS}$ of the neutron star, defined as the radius of a sphere with the same volume as the deformed rotating neutron star, "effective'' radius defined as $R_{\rm eff,NS} = \sqrt {I_{\rm NS}/M_{\rm NS}}$, and neutron star radii at the pole ( $R_{\rm NS,pol}$) and equator ( $R_{\rm NS,eq}$), for Model s15_64_r. Also shown is the eccentricity $\eta _{\rm NS} \equiv \sqrt {1-R_{\rm NS,pol}^2 / R_{\rm NS,eq}^2}$.

Open with DEXTER

As already discussed in earlier publications (Janka & Keil 1998; Janka et al. 2001; Keil 1997) the differential rotation, which is accounted for by the first term in Eq. (9), tends to damp convective activity near the polar axis in the PNS, see Figs. 26e,f. The same effect is found for the hot-bubble layer, see Figs. 26c,d. At larger distances from the polar axis the distribution of the specific angular momentum is flatter and convection is only weakly affected by rotation, i.e. convection develops in regions inside the PNS and between the gain radius and the shock basically similar to the situation in Model s15_32. PNS convection therefore has effects analogous to those discussed in Sect. 3.1. The damping of convection near the poles leads to a slightly slower effective transport of lepton number (and energy), see Fig. 28 for the situation after 200 ms post-bounce evolution, especially at later times when the convection near the axis is more strongly suppressed due to a steeper gradient of j_z, see Figs. 26e,f. Another new feature in the rotating model is the transport and redistribution of angular momentum by convection. Figure 25b shows how convection in the PNS produces a flat radial profile of j_z near the equator. This happens despite the stabilizing effect of the initially positive derivative of j_z because the negative entropy gradient dominates the Solberg-Høiland criterion and drives convective instability.

How does rotation affect the possibility for getting explosions by the delayed neutrino-heating mechanism? The analysis by Yamasaki & Yamada (2005) suggests that rotation can appreciably improve the conditions for shock revival in case the initial rotation frequency is at least $f \sim 0.1~$Hz at 1000 km. This corresponds to a rotation period of 10 s and is therefore only slightly faster than the rotation considered in Model s15_64_r (where the initial spin period is 12 s in the iron core). Yamasaki & Yamada (2005) found a sizable reduction by 25% of the "critical neutrino luminosity'' for starting a neutrino-driven explosion at a mass accretion rate of about 1$~M_\odot~$s^-1 or lower. The effects of rotation are, however, diverse and modify the structure of the collapsing star, the convection in the core, the gas motion behind the shock, and the radius and neutrino emission of the forming neutron star. A separate discussion of selected effects as done by Yamasaki & Yamada (2005) can therefore be misleading. In our simulations all effects of rotation are fully coupled and we can assess the question how these effects in combination determine the conditions for the delayed explosion mechanism. To this end we again compare the results of our rotating model with the non-rotating counterparts.

Rotation turns out to be helpful, but to a much lesser extent than estimated by Yamasaki & Yamada (2005). Certainly, the shock radius is significantly larger than in the non-rotating models of the 15$~M_\odot$ star (Fig. 15), the effective advection timescale  $\tau_{\rm adv}^\ast$ through the gain layer correspondingly becomes longer by up to a factor of about 4, and the mass in the gain layer, $M_{\rm gl}$, increases to a value that is two or three times larger than in the non-rotating models. These differences are on the one hand an indirect consequence of the structural changes of the PNS (cf. Sect. 3.1), on the other hand they result directly from the influence of centrifugal forces on the fluid flow in the gain layer. In contrast, however, the total heating rate, $\dot E_{{\rm gl}}$, and the ratio of advection to heating timescale increase only slightly (Fig. 17). The latter remains significantly below unity, because the heating timescale  $\tau_{\rm heat}$ has increased considerably due to the lower energies and luminosities of the neutrinos radiated from the rotating PNS (Figs. 12, 13). The rotating model therefore shows no tendency to develop an explosion until we stopped our simulation at about 280 ms after bounce.

If angular momentum conservation holds during the subsequent evolution, one can use our results at the end of the simulated evolution to roughly estimate the final angular frequency of the cold neutron star after its neutrino cooling and contraction phase, $\overline{\Omega}_{\rm cNS}$, assuming it to be a rigid rotator. The "NS'' is defined here as the mass at densities above $10^{11}~\ensuremath{{\rm g}~{\rm cm}^{-3}}$. In Fig. 29 we show the time evolution of several quantities for this such defined "NS'', namely of the radius $R_{\rm NS}$, defined as the radius of a sphere with the same volume as the deformed proto-neutron star, and of the average angular frequency, defined as $\overline {\Omega }_{\rm NS} = I_{\rm NS}/J_{\rm NS}$, where $I_{\rm NS}$ is the moment of inertia and $J_{\rm NS}$ is the total angular momentum of the "NS''. Taking the conditions at t=270 ms after bounce and assuming angular momentum convervation - no processes happen which transport angular momentum out of the NS or add mass and angular momentum to it -, the angular frequency of the NS after self-similar contraction (i.e., the shape of the star does not change and thus its moment of inertia scales with the square of the average NS radius) to a final radius $R_{\rm cNS}=10{-}12$ km is

 

$$\overline{\Omega}_{\rm cNS} \overline{\Omega}_{\rm NS,sim} \left(\frac{R_{\rm NS,sim}}{R_{\rm cNS}}\right)^2$$ =

$$\simeq 0.93 {\rm ~rad~ms}^{-1} \left(\frac{27 {\rm ~km}}{10\mbox{--}12{\rm ~km}}\right)^2 \simeq 5~...~7 {\rm ~rad~ms}^{-1},$$ (11)

where $\overline{\Omega}_{\rm NS,sim}$ and $R_{\rm NS,sim}$ are the average angular frequency and average radius of the "NS'' at the end of the simulation, respectively (Fig. 29). This corresponds to a spin period of $P_{\rm cNS}\equiv2\pi / \Omega_{\rm cNS} \simeq 1$ ms, which is at least 10-100 times faster than the typical birth period of pulsars estimated from observations (see Heger et al. 2005, and references therein). Therefore our Model s15_64_r is on the extreme side of conditions which are present in collapsing stellar cores in typical cases (see also the discussion in Ott et al. 2006).

   
3.5 Anisotropy of the neutrino emission

In this section we shall discuss the variations of the neutrino emission with polar angle, which are a consequence of PNS convection, convective overturn in the layer between gain radius and shock, and of rotation. For this purpose we will mostly concentrate on the non-rotating Model s112_128_f, which develops extremely strong convection and anisotropy, and on our case with rotation, Model s15_64_r. The anisotropy of the neutrino emission from rotating collapsing stellar cores was first evaluated by Janka & Mönchmeyer (1989b,a) analytically and by means of post-processing Monte Carlo transport calculations using 2D core-collapse models. They found that the neutrino emission is higher in the polar direction than in the equatorial plane, verifying von Zeipel's law of gravity darkening in case of neutrinos from supernova cores. This result was confirmed by Kotake et al. (2003), who also performed a post-processing analysis of 2D core-collapse calculations in which a trapping scheme for the neutrino treatment was employed. It must be pointed out, however, that none of these approaches was self-consistent: the emission anisotropy and the neutrino treatment in the hydrodynamic simulations were calculated with different approximative methods, the feedback of neutrino transport on the hydrodynamics was not taken into account (a trapping scheme only computes local source terms), and the evaluation was based on simplifying concepts, e.g., the assumption that a well-defined average neutrinosphere exists. The calculated anisotropies must therefore be interpreted with great caution and conclusions drawn on their basis may be very misleading, in particular concerning their implications for the explosion mechanism. The latter depends on many competing effects and a discussion requires a coupled and consistent treatment of neutrinos and hydrodynamics. Walder et al. (2005) recently performed such simulations of rotational core collapse, in which 2D flux-limited neutrino diffusion was followed self-consistently in 2D hydrodynamic models. These authors pointed out the fact that the neutrino emission is always more isotropic than the rotationally deformed or convectively distorted mass distribution because of the non-locality of the transport, i.e., the neutrino distribution at every point is a superposition of the irradiation from different contributing directions. The paper by Walder et al. (2005), however, does not provide quantitative information in a form which would allow us to make comparisons with their results. In the following we will therefore concentrate on the presentation and discussion of our models.

Figure 30: a) Lateral maxima of the $\nu _{{\rm e}}$ flux relative to the average $\nu _{{\rm e}}$ flux versus time for Models s15_64_r and s112_128_f, evaluated at 400 km radius for an observer at rest. b) Time integral of the "isotropic equivalent luminosity'' $L_{\rm iso}$ - calculated by assuming that the flux in one polar direction is representative for all other directions - versus $\cos\vartheta$ for different models and neutrinos. The integral was performed from 10 ms to 200 ms after bounce; $\cos\vartheta=\pm1$ correspond to the poles, $\cos\vartheta=0$ to the equator. The thin lines are for an evaluation at 400 km radius for an observer at rest, the thick lines are given only for $\nu _{{\rm e}}$ and show the result at the $\nu _{{\rm e}}$-sphere for an observer at rest.

Open with DEXTER

Figure 31: a) Ratio of electron neutrino flux to average flux (for an observer at rest at 400 km) versus cosine of the polar angle for Model s112_128_f for different post-bounce times. The times are picked such that large maxima of the flux ratio occur (see Fig. 30a). b) Same as panel a, but at the $\nu _{{\rm e}}$ sphere. The vertical solid lines mark the lateral boundaries of the region where the results were evaluated for Fig. 30a. Higher and lower latitudes were excluded because of possible artifacts from the reflecting boundaries.

Open with DEXTER

	Figure 32: Same as Fig. 31, but for Model s15_64_r. The pole of the rotating model corresponds to $\cos\vartheta = 1$.
Open with DEXTER

Figure 30a shows the lateral maximum of the ratio of the local neutrino flux to the average flux for  $\nu _{{\rm e}}$ as function of time. In case of Model s112_128_f spikes are visible with a typical width of about 5 ms and an amplitude that grows with time and reaches maximum values up to 2.4 towards the end of our simulation. These spikes are caused by "outbursts'' of neutrino radiation, which are associated with very narrow downflows of accreted matter, which reach deep into the cooling layer and create there local "hot spots'' (Fig. 21). These downflows and their location, however, are extremely time-dependent and nonstationary. This can be seen in Fig. 31a where the lateral flux variations are shown for three different times after core bounce. The large-amplitude flux variations are mostly a phenomenon which affects the neutrino emission from the cooling layer between gain radius and neutrinosphere. The amplitudes of $F_\vartheta/\left<F\right>$ are much smaller at the neutrinosphere (compare Figs. 31a and b).

Because the downflows are transient in space and time, the lateral variations disappear in the time-integrated flux at different latitudes, and only minor fluctations remain (Fig. 30b). Local maxima or minima at the lateral edges are a consequence of the reflecting boundary conditions, which cannot be penetrated but redirect the flow in the inward or outward direction.

The flux ratio maxima for Model s15_64_r also exhibit short-time variations with typical fluctuation timescales of about 5 ms, but with lower amplitudes (up to values of 1.5) than those in Model s112_128_f. This is mainly due to the fact that the convective layer in Model s15_64_r has a smaller radial extension, the downflows are less narrow and hit the cooling layer with less violence, producing less extreme local emission than in Model s112_128_f. An exception to this is the polar region of the rotating model, where the accretion flow is able to create a high-luminosity spot even at the neutrinosphere, causing larger latitudinal flux variations there than farther out (Figs. 32a,b).

While the time-integrated flux at large distances reveals a local maximum near the equator in Model s15_64_r (Fig. 30b), which is associated with a persistent equatorial downflow, the time-integrated flux at the neutrinosphere is nearly featureless with only a very shallow global pole-to-equator gradient. The latter is a consequence of the rotational flattening of the PNS, which, however, is too low to have significant effects on the instantaneous (Figs. 32b) or time-integrated (Fig. 30b) neutrinospheric emission.

We point out that a discussion of the lateral variation of the neutrino emission of our models is handicapped by the fact that our approximation of 2D neutrino transport disregards the lateral component of the neutrino flux vector and therefore tends to overestimate the angular asymmetry of neutrinos streaming out from radiating regions (for a more detailed discussion, see Buras et al. 2006). Truely multi-dimensional transport should therefore not only reveal smaller lateral variations of the time-integrated flux as our models do, but will also show less extreme angular variations of the instantaneous emission on short spatial wavelengths.

   
4 Summary and conclusions

We have presented results of a series of core-collapse and post-bounce simulations for different progenitor stars between 11.2$~M_\odot$ and 25$~M_\odot$, comparing 2D (axially symmetric) with 1D (spherically symmetric) calculations. Doing so, our main goals were (i) investigating the differences between convection in progenitors with different masses; (ii) exploring the effects of convection below the neutrinosphere ("PNS convection'') on the proto-neutron star structure, its neutrino emission, and the neutrino heating-layer behind the shock; (iii) investigating the role of hydrodynamic instabilities that affect the stalled accretion shock, i.e. convective overturn in the neutrino-heated "hot bubble'' layer ("HB convection'') and global low-mode nonradial instability of the accretion shock (termed SASI by Blondin et al. 2003; and possibly caused by the action of an advective-acoustic cycle according to Foglizzo 2001,2002); (iv) studying the effects of rotation; and (v) testing the influence of numerical aspects like the grid resolution, size of the angular wedge, and magnitude of seed perturbations for convection. Since our 2D neutrino-hydrodynamics code is a direct descendant of our 1D PROMETHEUS/VERTEX code, it is particularly well suited for performing such comparisons of 1D and 2D supernova models.

Convection inside the proto-neutron star starts 30-40 ms after bounce in all of our 2D models and encompasses a layer growing in mass until the end of our simulations (which were typically terminated about 250 ms after bounce). It leads to a more extended neutron star than in the 1D simulations with lower temperatures at the neutrinosphere. For this reason the mean energies of the neutrinos emitted from the neutrinosphere are reduced (up to 10% after 200 ms of PNS convection). Despite the larger radiating surface, the lower neutrinospheric temperatures also cause a slight reduction of neutrino luminosities during the first 150 ms after bounce. This holds in particular for  $\bar\nu_{{\rm e}}$, because the convective transport of lepton number maintains a higher electron degeneracy in the neutrinospheric region and accelerates the lepton number loss compared to 1D simulations. Only at $t\ga 150~$ms after bounce, convectively enhanced energy transport in the nascent neutron star also leads to increased energy loss, and the luminosities of heavy-lepton neutrinos become significantly (15%-20%) higher than in the spherical models.

PNS convection of the kind found in our simulations leads to a slightly reduced total energy deposition in the gain layer mainly because of the lower average energies of the radiated  $\nu _{{\rm e}}$ and  $\bar\nu_{{\rm e}}$. Since the effects of convection below the neutrinosphere are hard to disentangle from those of hydrodynamic instabilities in the neutrino-heating layer behind the stalled shock, we developed a simple "mixing algorithm''. It allowed us to reproduce all major effects of PNS convection in 1D simulations and thus to separate them from the consequences of multi-dimensional fluid flow in the postshock layer and to arrive at the above conclusion.

Convective overturn in the neutrino-heating layer remained rather weak in case of the 15$~M_\odot$ and 20$~M_\odot$ progenitors. The main reason for that is the rapid contraction of the accretion shock after its maximum expansion. This causes the gain layer to be very narrow and the infall velocities of the gas ahead and behind the shock to be very high. As a consequence, the advection timescale of the gas through the gain layer is very short compared to the typical neutrino-heating timescale. Buoyancy forces therefore hardly achieve bubble rise in the flow of gas accreted from the shock to the gain radius. As suggested by Janka & Keil (1998) and Janka et al. (2001) and verified by Thompson et al. (2005), the ratio of the advection timescale to the neutrino-heating timescale, $\tau _{\rm adv} / \tau _{\rm heat}$, turned out to be a useful diagnostic parameter to measure the proximity of a model to a neutrino-driven explosion. A necessary condition for an explosion is that the timescale ratio rises above unity for a time interval of at least the neutrino-heating timescale. In case of the 15$~M_\odot$ and 20$~M_\odot$ models, HB convection increases the heating rate and the timescale ratio to values only slightly larger than in the 1D simulations, but still roughly a factor of two below the critical limit. We therefore found explosions of these stars neither in spherical symmetry (in agreement with Liebendörfer et al. 2002; Thompson et al. 2003; Sumiyoshi et al. 2005) nor in 2D.

Also rotation did not change this negative outcome. We studied one 15$~M_\odot$ model with pre-collapse rigid iron core rotation of $\sim$12 s period, which leads to a neutron star with a spin period of about 1 ms, if the angular momentum of the core after collapse is conserved. Rotation of this size is probably on the extreme side of what can be expected for the cores of "normal'' supernovae, which are supposed to give birth to neutron stars with an initial period of 10 ms or more (see the discussions in Heger et al. 2005; and Ott et al. 2006). Our simulations reveal a number of important differences of the rotating model compared to its non-rotating counterparts. The proto-neutron star develops an eccentricity of more than 0.6 until we stopped the simulation at nearly 300 ms after bounce. At this time the luminosities of the radiated neutrinos are significantly smaller (10-20%) and their mean energies up to 2 MeV lower than in the non-rotating 2D model, because the equatorially more extended neutrinosphere is significantly cooler and energy is stored in rotation instead of being released by neutrinos. Despite the clear oblateness of the proto-neutron star, its rotation-induced emission anisotropy is very small. Nevertheless, rotation has a favorable influence on the conditions and parameters which determine neutrino-driven explosions. Centrifugal forces stabilize the accretion shock at larger radii, increase the advection timescale of the postshock gas significantly, and thus allow for a layer of well developed, strong convective overturn activity behind the shock. Because more mass stays in the gain layer for a longer time, the total energy deposition rate behind the shock is higher at later post-bounce times ( $t \ga 130~$ms p.b.) than in the non-rotating models. In spite of these healthy effects, however, the timescale ratio $\tau _{\rm adv} / \tau _{\rm heat}$remains still well below unity ( $\tau_{{\rm adv}}/\tau_{{\rm heat}} \la 0.6$).

Even without rotation postshock convection becomes violent in case of the 11.2 $~M_\odot$ star. The shock in this model is able to stay longer at large radii than in the more massive stars. This is due to the fact that the rather low-mass progenitor has a steeper density decline at the transition to the Si+O layer, which leads to a rapid decrease of the mass accretion rate of the shock at about 90 ms after bounce. This allows the shock to reexpand in adjusting to the situation of reduced ram pressure. The increased advection timescale gives convection the possibility to gain strength and thus to support the shock at a much larger radius than in the corresponding 1D model. Also the total neutrino heating rate behind the shock and the efficiency of net neutrino-energy transfer to the gas in the gain layer is higher by up to a factor of two. The timescale ratio $\tau _{\rm adv} / \tau _{\rm heat}$approaches unity and remains close to - but slightly below - this threshold until the end of our simulations. The 11.2 $~M_\odot$ model computed with a 90 $\hbox{$^\circ$ }$ lateral wedge therefore lingers at the border to success.

Such a situation is extremely sensitive to relatively little changes. We saw this when we repeated the simulation with a full 180 $\hbox{$^\circ$ }$ grid instead of using the wedge around the equator. While PNS convection turned out not to depend on the wedge size, convective activity in the neutrino-heating layer can change significantly when the available degrees of freedom are not constrained by periodic boundary conditions of a 90 $\hbox{$^\circ$ }$ equatorial wedge and therefore low-mode deformation of dipolar (l = 1) and quadrupolar (l = 2) character is allowed for. Convection becomes sufficiently strong so that the accretion shock continues to expand. This ensures that the effective advection timescale does not decrease after it has reached its maximum. At $t \ga 140~$ms after bounce, the timescale ratio $\tau _{\rm adv} / \tau _{\rm heat}$ then becomes larger than unity, thus further improving the conditions for efficient energy deposition by neutrinos in the postshock layer. After about 180 ms of post-bounce evolution the total energy in the gain layer becomes positive and continues rising because the mass in the gain layer and the energy per nucleon grow. The model has passed the critical threshold and is on its way to explosion. A closer inspection of the involved energies shows that this explosion is powered by neutrino heating.

This qualitative difference of the outcome of 2D simulations with 90 $\hbox{$^\circ$ }$ and 180 $\hbox{$^\circ$ }$ grids is another confirmation of the proximity of our 2D simulations, and in particular of the 11.2 $~M_\odot$ case, to a success of the convectively supported neutrino-driven mechanism. Together with the recent models for stars in the $8{-}10~M_\odot$ range with O-Ne-Mg cores, which explode even in spherical symmetry (Kitaura et al. 2006), our current results seem to indicate that the neutrino-heating mechanism is viable at least for stars near the low-mass end of supernova progenitors.

The sensitivity to numerical variations, however, also stresses the need to remove some of the shortcomings and limitations of axially symmetric simulations. One must suspect that in 3D simulations morphological differences of the structures (plumes instead of azimuthal tori), different growth rates of instabilities, or additional degrees of freedom (e.g. triaxial asymmetries and vortex motion caused by Coriolis forces) might lead to sizable quantitative differences which could be crucial when collapsing stellar cores are close to the threshold for explosion. Also the existence of the polar axis of a spherical or cylindrical coordinate grid is a potential source of numerical uncertainties, because it is a coordinate singularity which is impenetrable for approaching fluid flow and thus defines a preferred grid direction.

Our results therefore suggest the need to strive for 3D simulations, preferentially without the disadvantages connected with the polar grid axis. The importance of low-mode convection or low-mode hydrodynamical instabilities as suggested by our results implies that such simulations will have to be done for the full star and cannot be contrained to a limited wedge.

Acknowledgements

We are grateful to Almudena Arcones, Francisco Kitaura, Andreas Marek, and Leonhard Scheck for many helpful discussions and to an anonymous referee for careful reading and a long list of useful comments. Support by the Sonderforschungsbereich 375 on "Astro-Particle Physics'' of the Deutsche Forschungsgemeinschaft is acknowledged. The computations were performed on the NEC SX-5/3C and the IBM p690 "Regatta'' system of the Rechenzentrum Garching, and on the Cray T90 and IBM p690 "Jump'' of the John von Neumann Institute for Computing (NIC) in Jülich.

References

• Arcones Segovia, A. 2003, Diploma thesis, Technische Universität München (In the text)
• Blondin, J. M., & Mezzacappa, A. 2006, ApJ, 642, 401 [NASA ADS] [CrossRef]
• Blondin, J. M., Mezzacappa, A., & DeMarino, C. 2003, ApJ, 584, 971 [NASA ADS] [CrossRef]
• Bruenn, S. W., & Dineva, T. 1996, ApJ, 458, L71 [NASA ADS] (In the text)
• Bruenn, S. W., & Mezzacappa, A. 1997, Phys. Rev. D, 56, 7529 [NASA ADS] [CrossRef] (In the text)
• Bruenn, S. W., Raley, E. A., & Mezzacappa, A. 2004, ApJ,submitted [arXiv:astro-ph/0404099] (In the text)
• Buras, R., Rampp, M., Janka, H.-T., & Kifonidis, K. 2003, Phys. Rev. Lett., 90, 241101 [NASA ADS] [CrossRef]
• Buras, R., Rampp, M., Janka, H.-T., & Kifonidis, K. 2006, A&A, 447, 1049 [EDP Sciences] [NASA ADS] [CrossRef] (Paper I)
• Burrows, A., Hayes, J., & Fryxell, B. A. 1995, ApJ, 450, 830 [NASA ADS] [CrossRef] (In the text)
• Burrows, A., Livne, E., Dessart, L., Ott, C. D., & Murphy, J. 2006, ApJ, 640, 878 [NASA ADS] [CrossRef] (In the text)
• Cardall, C. Y. 2004, in Numerical Methods for Multidimensional Radiative Transfer Problems (RadConf2003), Proc. Int. Conf. in Heidelberg, Germany, Sept. 24-26, 2003, ed. R. Rannacher, & R. Wehrse [arXiv:astro-ph/0404401] (In the text)
• Cardall, C. Y., Razoumov, A. O., Endeve, E., & Mezzacappa, A. 2005, in Open Issues in Understanding Core Collapse Supernovae, proc. conf. at National Institute for Nuclear Theory, Univ. Washington, June 22-24, 2005, World Scientific, in press [arXiv:astro-ph/0510704] (In the text)
• Chandrasekhar, S. 1961, Hydrodynamics and Hydromagnetic Instability (New York: Dover Publications), 231 (In the text)
• Dessart, L., Burrows, A., Livne, E., & Ott, C. D. 2005, ApJ, submitted [arXiv:astro-ph/0510229] (In the text)
• Foglizzo, T. 2001, A&A, 368, 311 [EDP Sciences] [NASA ADS] [CrossRef]
• Foglizzo, T. 2002, A&A, 392, 353 [EDP Sciences] [NASA ADS] [CrossRef]
• Foglizzo, T., Scheck, L., & Janka, H.-T. 2005, A&A, submitted [arXiv:astro-ph/0507636] (In the text)
• Fryer, C. L. 1999, ApJ, 522, 413 [NASA ADS] [CrossRef] (In the text)
• Fryer, C. L., & Heger, A. 2000, ApJ, 541, 1033 [NASA ADS] [CrossRef] (In the text)
• Fryer, C. L., & Warren, M. S. 2002, ApJ, 574, L65 [NASA ADS] [CrossRef]
• Fryer, C. L., & Warren, M. S. 2004, ApJ, 601, 391 [NASA ADS] [CrossRef]
• Fryer, C. L., Benz, W., & Herant, M. 1996, ApJ, 460, 801 [NASA ADS] [CrossRef] (In the text)
• Heger, A., Langanke, K., Martínez-Pinedo, G., & Woosley, S. E. 2001, Phys. Rev. Lett., 86, 1678 [NASA ADS] [CrossRef]
• Heger, A., Woosley, S. E., & Spruit, H. C. 2005, ApJ, 626, 350 [NASA ADS] [CrossRef] (In the text)
• Herant, M., Benz, W., Hix, W. R., Fryer, C. L., & Colgate, S. A. 1994, ApJ, 435, 339 [NASA ADS] [CrossRef] (In the text)
• Janka, H.-T. 2001, A&A, 368, 527 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Janka, H.-T., & Keil, W. 1998, in Supernovae and Cosmology, ed. L. Labhardt, B. Binggeli, & R. Buser, Astronomisches Institut, Universität Basel, 7 [arXiv:astro-ph/9709012] (In the text)
• Janka, H.-T., & Mönchmeyer, R. 1989a, A&A, 209, L5 [NASA ADS]
• Janka, H.-T., & Mönchmeyer, R. 1989b, A&A, 226, 69 [NASA ADS]
• Janka, H.-T., & Müller, E. 1996, A&A, 306, 167 [NASA ADS] (In the text)
• Janka, H.-T., Kifonidis, K., & Rampp, M. 2001, in Physics of Neutron Star Interiors, ed. D. Blaschke, N. K. Glendenning, & A. Sedrakian, Notes in Physics, 578 (Berlin: Springer), 333 [arXiv:astro-ph/0103015] (In the text)
• Janka, H.-T., Buras, R., Kitaura Joyanes, F. S., et al. 2005, Nuclear Physics A, 758, 19c [NASA ADS] [CrossRef] (In the text)
• Kachelrieß, M., Tomàs, R., Buras, R., et al. 2005, Phys. Rev. D, 71, 063003 [NASA ADS] [CrossRef] (In the text)
• Kane, J., Arnett, D., Remington, B. A., et al. 2000, ApJ, 528, 989 [NASA ADS] [CrossRef] (In the text)
• Keil, W. 1997, Ph.D. Thesis, Technische Universität München (In the text)
• Keil, W., & Janka, H.-T. 1995, A&A, 296, 145 [NASA ADS] (In the text)
• Keil, W., Janka, H.-T., & Müller, E. 1996, ApJ, 473, L111 [NASA ADS] [CrossRef] (In the text)
• Kitaura, F. S., Janka, H.-T., & Hillebrandt, W. 2006, A&A, 450, 345 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Kotake, K., Yamada, S., & Sato, K. 2003, ApJ, 595, 304 [NASA ADS] [CrossRef] (In the text)
• Kotake, K., Sawai, H., Yamada, S., & Sato, K. 2004, ApJ, 608, 391 [NASA ADS] [CrossRef]
• Lai, D., & Goldreich, P. 2000, ApJ, 535, 402 [NASA ADS] [CrossRef]
• Langanke, K., Martínez-Pinedo, G., Sampaio, J. M., et al. 2003, Phys. Rev. Lett., 90, 241102 [NASA ADS] [CrossRef]
• Lattimer, J., & Swesty, F. 1991, Nucl. Phys. A, 535, 331 [NASA ADS] [CrossRef] (In the text)
• Liebendörfer, M., Mezzacappa, A., Thielemann, F., et al. 2001, Phys. Rev. D, 63, 103004 [NASA ADS] [CrossRef]
• Liebendörfer, M., Messer, O. E. B., Mezzacappa, A., et al. 2002, in Proceedings of the 11th Workshop on Nuclear Astrophysics, Ringberg Castle, Tegernsee, Germany, February 11-16, 2002, ed. E. Müller, & W. Hillebrandt, 126 [arXiv:astro-ph/0203260] (In the text)
• Liebendörfer, M., Mezzacappa, A., Messer, O. E. B., et al. 2003, Nucl. Phys. A, 719, C144 [NASA ADS] [CrossRef] (In the text)
• Liebendörfer, M., Messer, O. E. B., Mezzacappa, A., et al. 2005, ApJS, 150, 263 [NASA ADS] [CrossRef]
• Limongi, M., Straniero, O., & Chieffi, A. 2000, ApJS, 129, 625 [NASA ADS] [CrossRef] (In the text)
• Livne, E., Burrows, A., Walder, R., Lichtenstadt, I., & Thompson, T. A. 2004, ApJ, 609, 277 [NASA ADS] [CrossRef] (In the text)
• Marek, A., Janka, H.-T., Buras, R., Liebendörfer, M., & Rampp, M. 2005a, A&A, 443, 201 [EDP Sciences] [NASA ADS] [CrossRef]
• Mezzacappa, A., & Blondin, J. M. 2003, in From Twilight to Highlight: The Physics of Supernovae, ed. W. Hillebrandt, & B. Leibundgut (Berlin: Springer), 63 (In the text)
• Mezzacappa, A., Calder, A. C., Bruenn, S. W., et al. 1998, ApJ, 495, 911 [NASA ADS] [CrossRef] (In the text)
• Müller, E., Rampp, M., Buras, R., & Janka, H.-T. 2004, ApJ, 603, 221 [NASA ADS] [CrossRef] (In the text)
• Nomoto, K., & Hashimoto, M. 1988, Phys. Rep., 163, 13 [NASA ADS] [CrossRef]
• Ohnishi, N., Kotake, K., & Yamada, S. 2006, ApJ, 641, 1018 [NASA ADS] [CrossRef]
• Ott, C. D., Burrows, A., Livne, E., & Walder, R. 2004, ApJ, 600, 834 [NASA ADS] [CrossRef]
• Ott, C. D., Burrows, A., Thompson, T. A., Livne, E., & Walder, R. 2006, ApJS, 164, 130 [NASA ADS] (In the text)
• Rampp, M., & Janka, H.-T. 2002, A&A, 396, 361 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Scheck, L. 2005, Ph.D. Thesis, Technische Universität München, in preparation (In the text)
• Scheck, L., Plewa, T., Janka, H.-T., Kifonidis, K., & Müller, E. 2004, Phys. Rev. Lett., 92, 011103 [NASA ADS] [CrossRef]
• Sumiyoshi, K., Yamada, S., Suzuki, H., et al. 2005, ApJ, 629, 922 [NASA ADS] [CrossRef] (In the text)
• Swesty, F. D., & Myra, E. S. 2005a, in Open Issues in Understanding Core Collapse Supernovae, Proc. Conf. at National Institute for Nuclear Theory, Univ. Washington, June 22-24, 2005 (World Scientific) [arXiv:astro-ph/0506178]
• Swesty, F. D., & Myra, E. S. 2005b, to appear in Journal of Physics: Conf. Ser.; talk at SciDAC 2005, San Francisco, June 26-30, 2005 [arXiv:astro-ph/0507294]
• Tassoul, J.-L. 1978, Theory of Rotating Stars (Princeton: Princeton University Press)
• Thompson, C. 2000, ApJ, 534, 915 [NASA ADS] [CrossRef]
• Thompson, T. A., Burrows, A., & Pinto, P. A. 2003, ApJ, 592, 434 [NASA ADS] [CrossRef] (In the text)
• Thompson, T. A., Quataert, E., & Burrows, A. 2005, ApJ, 620, 861 [NASA ADS] [CrossRef] (In the text)
• Walder, R., Burrows, A., Ott, C. D., et al. 2005, ApJ, 626, 317 [NASA ADS] [CrossRef] (In the text)
• Wilson, J. R., & Mayle, R. 1988, Phys. Rep., 163, 63 [NASA ADS] [CrossRef]
• Wilson, J. R., & Mayle, R. 1993, Phys. Rep., 227, 97 [NASA ADS] [CrossRef]
• Woosley, S. E., & Weaver, T. A. 1995, ApJS, 101, 181 [NASA ADS] [CrossRef]
• Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Rev. Mod. Phys., 74, 1015 [NASA ADS] [CrossRef] (In the text)
• Yamasaki, T., & Yamada, S. 2005, ApJ, 623, 1000 [NASA ADS] [CrossRef] (In the text)

 

  
Online Material

   
Appendix A: The progenitor models

Table A.1: List of progenitors used in the simulations. $M_{{\rm Fe-Si}}$ and $M_{{\rm Si-O}}$ are the enclosed masses at composition interfaces defined by entropy jumps. With (e) we denote that a shell interface is connected with a gradual enrichment of the lighter nucleus, i.e. these cases are Fe-FeSi and Si-SiO interfaces where the mass fractions of Si or O grows gradually outwards. The underlined numbers indicate an entropy increase of more than $1~\ensuremath{k_{\rm B}}$ per nucleon in case of an Fe-Si interface and an entropy increase of more than $2~\ensuremath{k_{\rm B}}$/by in case of a Si-O interface.

Figure A.1: The progenitor data for temperature T, electron fraction ${Y_{\rm e}}$, density $\rho$, and the mass fractions of Si and O, $X_{\rm Si}$ and $X_{\rm O}$, respectively, as provided to us by the stellar evolution modelers (see Table A.1 for references). The binding energy and entropy were derived with the EoS used in our simulations.

The properties of the nine progenitor models used in this work are summarized in Fig. A.1 and Table A.1. In Fig. A.1 the initial pre-collapse data as given by the stellar evolution modelers are displayed, in Figs. A.2 and A.3 we compare the cores at a central density of $\rho_{\rm c} = 10^{11}~\ensuremath{{\rm g}~{\rm cm}^{-3}}$, which was reached by evolving the models with our 1D Boltzmann transport code VERTEX. Up to this point the infall velocities are still subsonic and neutrinos stream off almost unhindered. Compared to the original progenitor data in Fig. A.1 the electron number ${Y_{\rm e}}$ has changed significantly and strongest for those models which started out with the lowest initial densities at the center. The evolution proceeded nearly adiabatically, the entropy changes are therefore small.

Looking at Fig. A.2, we see that the density structure of the inner $\sim$ $1.0~\ensuremath{M_{\odot}}$ of the iron core is extremely similar in all progenitors, correlated with only small differences in the electron fraction  ${Y_{\rm e}}$ (Fig. A.3). The central value of ${Y_{\rm e}}$ varies only by 5% between the progenitors. In contrast, the entropy per baryon s and infall velocity v exhibit differences of up to 40%. There is a general trend that |v| and s increase with the ZAMS mass, while  ${Y_{\rm e}}$ decreases.

Figure A.2: Radial profiles of the models during core collapse when the central density is $\rho_{\rm c} = 10^{11}~\ensuremath{{\rm g}~{\rm cm}^{-3}}$.

Figure A.3: Profiles as functions of enclosed mass of the models during core collapse when the central density is $\rho_{\rm c} = 10^{11}~\ensuremath{{\rm g}~{\rm cm}^{-3}}$. The data are only shown for the stellar mass range that was included in the simulations. It corresponds to $r \leq 10~000$ km.

Outside of the core, larger differences between the progenitors exist and are associated with the location of the interfaces between layers of different chemical composition and the density structure in these layers. In particular, the Fe-Si and the Si-O interfaces can have significant influence on the evolution of the supernova shock. The enclosed mass at which the interfaces are located differs strongly between the progenitors and increases non-monotonically with the ZAMS mass, see Table A.1. In most cases the composition changes discontinuously from the heavy to the lighter nucleus at the interfaces, but also a more gradual enrichment of a layer with elements of the neighboring shell is possible. In both cases the progenitor structure shows a more or less large entropy jump at the interfaces. For larger entropy jumps, which are underlined in Table A.1, also large steps occur in the density profile.

   
Appendix B: One-dimensional models: more details

Table B.1: Characteristic parameters of the 1D models for the phases of collapse, prompt shock propagation, and neutrino burst. $t_{\rm coll}$ is the time between the moment when the collapsing core reaches a central density of  $10^{11}~\ensuremath{{\rm g}~{\rm cm}^{-3}}$ and the moment of shock creation (which is nearly identical with the time of core bounce), i.e. the time when the entropy behind the shock first reaches a value of  $3~\ensuremath{k_{\rm B}}$/by. The shock creation radius  $r_{\rm sc}$ and enclosed mass  $M_{\rm sc}$ are defined by the radial position where this happens. The energy loss $E_{\rm coll}^{\nu ,{\rm loss}}$ via neutrinos during the collapse phase is evaluated by integrating the total neutrino luminosity (for an observer at rest at r=400 km) over time from the moment when the core reaches $\rho_{\rm c} = 10^{11}~\ensuremath{{\rm g}~{\rm cm}^{-3}}$ until the moment when the dip in the $\nu _{{\rm e}}$ luminosity is produced around 2.5 ms after shock formation. We call the time when the velocities behind the shock drop below $10^7~~\ensuremath{{\rm cm}~{\rm s}^{-1}}$ the end of the prompt shock propagation phase. This time, $t_{\rm pse}$, is measured relative to the moment of shock creation. At the end of the prompt shock propagation phase, the shock is at radius $r_{\rm pse}$ and its enclosed mass is $M_{\rm pse}$. The time of the $\nu _{{\rm e}}$ burst, $t_{{\nu _{\rm e}}-{\rm burst}}$, is defined as the post-bounce time when the $\nu _{{\rm e}}$ luminosity maximum is produced at the shock, which then is located at the radius $r_{{\nu _{\rm e}}-{\rm burst}}$ and mass $M_{{\nu _{\rm e}}-{\rm burst}}$. Finally, the energy emitted during the prompt $\nu _{{\rm e}}$ burst is defined as the time-integrated $\nu _{{\rm e}}$ luminosity for the FWHM of the burst, $E_{\rm burst}^\nu$, evaluated at 400 km for an observer at rest.

With the exception of Model n13, the collapse times from a core central density of $\rho_{\rm c} = 10^{11}~\ensuremath{{\rm g}~{\rm cm}^{-3}}$ until bounce differ only by 10-20% between the models, lying between $\sim$20 and $\sim$25 ms (see Fig. B.1 and Table B.1). In general, the results listed in Table B.1 and visible in Figs. A.3 and B.1-B.3 reveal an astonishing degree of convergence of the core evolution for the different progenitors during the phases of collapse, shock formation, prompt shock propagation and breakout $\nu _{{\rm e}}$burst (cf. also Liebendörfer et al. 2002). Initially, the stellar cores differ in the amount of deleptonization and have been evolved differently close to the onset of collapse. Less deleptonized cores, however, need a longer time until the collapse becomes dynamical (because their ${Y_{\rm e}}$ and thus electron pressure are higher and photodisintegrations of Fe-group nuclei proceed more slowly), thus allowing deleptonization to catch up with that of the more evolved progenitors. Despite of remaining differences at $\rho_{\rm c} = 10^{11}~\ensuremath{{\rm g}~{\rm cm}^{-3}}$ (Figs. A.3 and B.2), the central quantities and radial profiles, in particular of entropy and ${Y_{\rm e}}$, become very similar after neutrino trapping sets in at $\rho_{\rm c} \ga 10^{12}~\ensuremath{{\rm g}~{\rm cm}^{-3}}$(Figs. B.2, B.3). This suggests a strongly self-regulating character of the hydrodynamics coupled with the neutrino transport, which ensures that the characteristic properties of the homologous core at bounce are almost independent of the initial conditions (Liebendörfer et al. 2002). Only Model s25a28 has such a high value of the entropy and such a low value of  ${Y_{\rm e}}$ at the time when $\rho_{\rm c} = 10^{11}~\ensuremath{{\rm g}~{\rm cm}^{-3}}$ (Figs. A.3 and B.2) that both quantities do not fully converge to those of the other models. In this model the higher entropy implies that the EoS yields a much larger free proton abundance, which leads to slightly stronger deleptonization and thus a somewhat larger entropy (the entropy level is higher by about $\Delta s\sim 0.07~\ensuremath{k_{\rm B}} /{\rm nucleon}$) after trapping.

As a consequence of the similar collapse evolution, all models lose approximately the same amount of energy via neutrino emission during collapse, i.e. around 10⁵¹ erg, and because of the similar core structure and collapse history, the shock in all cases is created at an enclosed mass of $M_{\rm sc} \simeq 0.49\pm 0.01~\ensuremath{M_{\odot}}$, corresponding to a radius of $r_{\rm sc} \simeq 10.7\pm 0.1~{\rm km}$ (Table B.1 and Fig. B.3a). Following Bruenn & Mezzacappa (1997) we define the shock creation (sc) time and location by the moment and position where the entropy first reaches $3~\ensuremath{k_{\rm B}}$ per baryon. From now on all times will be normalized to the time of core bounce, $t_{\rm cb}$.

Also the prompt shock propagation is quite similar. In most models, the velocities behind the shock become negative after 0.9 ms (time  $t_{\rm pse}$ in Table B.1) when the shock encloses a mass of  $M_{\rm pse}=0.78\pm 0.01~\ensuremath{M_{\odot}}$ and has reached a radius of  $r_{\rm pse} = 32\pm 1{\rm ~km}$ (pse stands for "prompt shock ends''). Only in case of the progenitors from Limongi et al. (2000) the prompt shock pushes out a bit farther. It is interesting to note that in all cases the stagnation of the shock happens earlier than the  $\nu _{{\rm e}}$ burst is released. The corresponding time $t_{{\nu _{\rm e}}-{\rm burst}}$ in Table B.1 is defined as the moment when the $\nu _{{\rm e}}$ luminosity maximum is produced at the shock.

Although the prompt shock "stalls'' in the above defined sense, the radial expansion velocity of the shock remains large, see Fig. B.4, upper left panel, because of mass of the collapsing star being accreted through the shock and accumulating on the central, collapsed core. When the shock has passed the neutrinosphere after 4 ms at a radius of 60-70 km (Table B.1), the energy and lepton number drain via neutrino emission in the $\nu _{{\rm e}}$ burst reduces the thermal and degeneracy pressure of the electrons in the accreted material, so that the shock continues to expand more slowly. The strength of this neutrino burst is very similar for all models, see Fig. B.4, lower left panel (cf. also Kachelrieß et al. 2005).

Our models show that the four parameters which characterize the subsequent quasi-stationary accretion phase, the mass accretion rate through the shock, $\partial_t M_{\rm sh}$, the mass $M_{\rm PNS}$ and radius $r_{\rm PNS}$ of the proto-neutron star (PNS), and the neutrino luminosity $L_\nu$, are not independent but coupled (Fig. B.4). The governing variable is the time-dependent mass accretion rate through the shock, which is determined by the progenitor structure. Since the mass accreted by the shock is further advected onto the PNS with a small time delay, the mass of the PNS is essentially the time integral of $\partial_t M_{\rm sh}$ plus an initial value. Starting at $t_{\rm pb} = 4{\rm ~ms}$, this initial PNS mass is approximately $1.0~\ensuremath{M_{\odot}}$ for all models due to the similar progenitor core structure. Also the neutrino luminosity depends on  $\partial_t M_{\rm sh}$: the gravitational binding energy of the accreted matter must be radiated away when accretion proceeds in a stationary way. The total neutrino luminosity therefore contains a part from the cooling of the contracting core of the PNS (which loses its binding energy over a timescale of about 10 s) plus a contribution from the settling accretion layer, which is proportional to  $\partial_t M_{\rm sh}$ with some lag, because the material has to fall from the shock to the cooling layer, where neutrinos are then released over a thermal cooling timescale of some 10 ms. Interestingly, also the neutrinosphere radii reveal a variability with  $\partial_t M_{\rm sh}$. A large (small) accretion rate leads to an approximately stationary solution with more (less) matter piling up on the nascent neutron star before it can radiate away its energy in form of neutrinos. Thus the PNS obtains a hot, extended (cool, narrow) mantle with relatively high (low) densities and optical depths, releasing neutrinos only on a longer (shorter) timescale. The neutrinosphere therefore moves to a larger (smaller) radius $r_\nu$ in case of mass accretion proceeding at high (low) rates (Fig. B.4).

Figure B.1: Central density versus time remaining until core bounce for all 1D models.

Figure B.2: Central entropy (panel  a)) and central electron and lepton fraction (panel  b)) versus central density during core collapse for all 1D models. The vertical dotted line marks the time of comparison at a central density  $10^{11}~\ensuremath{{\rm g}~{\rm cm}^{-3}}$ in Figs. A.2 and A.3.

Figure B.3: Entropy (panel  a)) and electron and lepton fraction (panel  b)) versus enclosed mass at the time of shock formation.

In summary, all the variables which determine the structure of the postshock accretion layer depend on the time evolution of  $\partial_t M_{\rm sh}$: a high $\partial_t M_{\rm sh}$ leads to larger values of $r_\nu$ and $L_\nu$, and a faster increase of  $M_{\rm PNS}$. Figure B.4 shows this dependence. Looking, e.g., at $t_{\rm pb} = 80{\rm ~ms}$ we see that $\partial_t M_{\rm sh}$ varies by a factor of five between the models, being 0.8 and 1.3  $\ensuremath{M_{\odot}} ~{\rm s}^{-1}$ in case of Models n13 and s11.2, respectively, at the low mass end of considered progenitors, and 3-4  $\ensuremath{M_{\odot}} /{\rm s}$ for Models l25 and s25a28 at the high end of the progenitor mass spectrum. Consistently, $r_\nu$ at $t_{\rm pb} = 80{\rm ~ms}$ is only 60 km for the low-mass progenitors and 25% higher, i.e. about 80 km, for the high-mass models. Also the luminosities differ by a factor of two, and the PNS masses show differences of order 25%.

The shock radius follows the time evolution of the neutron star radius qualitatively, but with some time delay. Interestingly, during phases of slowly changing  $\partial_t M_{\rm sh}$ the influence of different values of  $\partial_t M_{\rm sh}$ on the shock position is rather modest. The shock radii are nearly identical for all progenitors until $\sim$25 ms after bounce and differ at most by about 10% until about 80 ms post-bounce. During this phase, the significantly larger differences in  $\partial_t M_{\rm sh}$, $M_{\rm PNS}$, $r_\nu$, and $L_\nu$ mentioned above seem to partly compensate each other in their influence on the shock radius. For higher values of  $\partial_t M_{\rm sh}$ the shock radius is smaller because of the higher ram pressure. A larger  $M_{\rm PNS}$ increases the gravitational pull and thus also lowers  $r_{\rm sh}$. But on the other hand a higher $L_\nu$ increases the heating and thus leads to an expansion of the heating region, which is supported from below by a larger PNS (because $r_\nu$ is larger).

Of course, the counteracting effects do not compensate each other perfectly. As can be seen at later times, the shock radii differ more strongly when the PNSs begin to show a wider spread in mass, although with nearly the same radius, and in particular when the mass accretion rates show differences of a factor of about 10 instead of the factor of 2 at early times. The analytic study of Arcones Segovia (2003) and Arcones & Janka (in preparation) reveals that the steady-state accretion shock radius is a sensitive function of the PNS radius ( $r_{\rm sh} \propto r_\nu^{1.5\dots2}$), and decreases less strongly with higher mass accretion rate and proto-neutron star mass (see also Fryer et al. 1996).

Figure B.4: Post-bounce evolution of the models in spherical symmetry. The left panels show phases around the core bounce with the last stage of infall (t < 0), prompt shock propagation, neutrino burst, and early accretion phase of the shock. The right panels depict the post-bounce accretion phase of the stalled shock. From top to bottom, the panels display the shock radius, the mass accretion rate through the shock, the mass enclosed by the shock ( left) or PNS baryonic mass ( right), the PNS radius, and the electron neutrino ( left) or total neutrino ( right) luminosity at r=400 km for an observer at rest. The PNS mass and radius are defined by the corresponding values at the $\nu _{{\rm e}}$ sphere.

Note that a sudden decrease of  $\partial_t M_{\rm sh}$ as it happens when an infalling composition shell interfaces reaches the shock (see Sect. 2, also nicely demonstrates the dependence of $r_\nu$, $L_\nu$, and  $M_{\rm PNS}$ on this variable: immediately after such a drop, $L_\nu$ starts decreasing, the decrease of $r_\nu$ becomes steeper, and the increase of  $M_{\rm PNS}$ slows down (Fig. B.4). Also, the total lepton number and energy loss via neutrino emission, see Fig. B.5, depends strongly on  $\partial_t M_{\rm sh}$ and thus on the progenitor structure during the post-bounce accretion phase. These facts could be used to infer the progenitor structure from a supernova neutrino signal which is measured with good time resolution (cf. also Liebendörfer et al. 2003).

Figure B.5: Integrated lepton number ( top) and energy loss of the 1D models versus time.

   
Appendix C: Linear analysis of PNS structural changes

A linear analysis is performed to determine the sensitivity of the PNS structure, i.e. contraction or expansion, to changes or redistribution of lepton number and energy in some layer. This analysis is then applied to the main effects of convective activity in the nascent neutron star, (1) the transport of lepton number and energy within the convective layer, and (2) the differences in the loss of lepton number and energy by neutrino emission in 2D simulations compared to 1D models.

We demonstrate the analysis in case of lepton number variations, holding the specific internal energy fixed, but variations of the internal energy can be investigated in the same way. Assuming Newtonian gravity, which is sufficient here because we are interested in qualitative, not quantitative results, the hydrostatic structure of the PNS is determined by

$$\frac{{\rm d}p(r)~}{{\rm d}r~} = - \frac{GM(r)\rho}{r^2} ,$$ (C.1)

with G being the Newtonian gravitational constant. Using the relation ${\rm d}M~ = 4\pi\rho r^2{\rm d}r~$, we get

$$\frac{{\rm d}p(M)~}{{\rm d}M~} = - \frac{GM}{4\pi r^4(M)}\cdot$$ (C.2)

An infinitesimal change $\delta Y(M)$ of lepton number, e.g. by removing one electron from a layer with infinitesimal mass $\Delta M$and baryon number  $\Delta N_{\rm by}$ at mass coordinate $M_{\rm p}$ (we here and in the following drop the subscript of  ${Y_{\rm lep}}$), i.e.

$$\delta Y(M) = -1\cdot \Delta N_{\rm by}^{-1} = -\left(\Delta M/m_{\rm by}\right)^{-1}$$ (C.3)

if $M_{\rm p} < M < M_{\rm p} + \Delta M \nonumber$, and $\delta Y(M) = 0$ otherwise, will lead to a variation of the pressure according to

$$\delta p(M) = \left.\frac{\partial p}{\partial Y}\right\vert ... ...artial p}{\partial \rho}\right\vert _Y (M) ~\delta \rho (M) ,$$ (C.4)

where $\delta \rho (M)$ is the structural change induced by the change of lepton number. The new hydrostatic solution is

$$\frac{{\rm d}(p+\delta p)~}{{\rm d}M~} = - \frac{GM}{4\pi (r+\delta r)^4} ,$$ (C.5)

where $\delta r(M)$ is the movement of mass shell M.

Expanding the RHS to first order, inserting Eq. (C.4), and subtracting Eq. (C.2) leads to

$$\frac{{\rm d}~}{{\rm d}M~} \left[ \left.\frac{\partial p}{\p... ... \cdot \delta \rho \right] = \frac{GM~\delta r}{\pi~r^5} \cdot$$ (C.6)

Now setting the RHS of Eq. (C.6) to zero for the moment, the solution reads

$$\delta \rho(M) = - \left.\frac{\partial p}{\partial Y}\righ... ...tial \rho}\right\vert _Y (M) \right]^{-1} \cdot \delta Y(M) .$$ (C.7)

In terms of the "mass shift'' $\delta M$ the resulting expansion/contraction of the PNS is now evaluated at a radius $R > r(M_{\rm p})$:

$$\delta M(R) \equiv \int_0^R 4\pi \delta \rho(M(r')) r'^2 {\rm d}r'~$$

$$\simeq \left.\frac{\partial p}{\partial Y}\right\vert _\rho (M_{\rm p}) ... ...ght\vert _Y (M_{\rm p}) \right]^{-1} \cdot \frac{m_{\rm by}}{\rho(M_{\rm p})} ,$$ (C.8)

where we have expanded to first order and made use of Eqs. (C.7) and (C.3) and of the relation  $\Delta M = 4\pi r^2\rho(M_{\rm p})\delta r$. The change in mass $\delta M$is independent of R and means expansion, if negative (i.e. less mass within a given radius R).

A relation between changes of mass and radius can be obtained by using $M = \int_0^{r} 4\pi \rho r'^2 {\rm d}r'~ \equiv \int_0^{r+\delta r} 4\pi (\rho + \delta \rho) r'^2 {\rm d}r'~$. Expanding to first order, we thus get a relation between $\delta r$ and  $\delta \rho$:

$$\delta r = - \frac{1}{4\pi\rho r^2} \int_0^r 4\pi\delta\rho r'^2 {\rm d}r'~ .$$ (C.9)

Note that the radius displacement $\delta R(M)$ depends on R, see Eq. (C.9), and is therefore less suitable for our analysis.

In order to take into account the RHS of Eq. (C.6) we apply an iterative procedure: The solution $\delta r_i(M)$ of iteration step i is inserted on the RHS of Eq. (C.6), which is then numerically integrated over M to find the solution $\delta \rho_{i+1} (M)$. In the first iteration step the solution for $\delta \rho (M)$ from Eq. (C.7) is used to calculate $\delta r_i(M)$ from Eq. (C.9). We find that for typical PNS profiles Eq. (C.7) provides a very good approximation already; the next iteration step changes $\delta M$ by 10^-5 at most.

Figure C.1: Results of a linear analysis of the structural changes of the PNS in response to the removal of lepton number or energy, respectively, at different densities and thus enclosed masses $M_{\rm p}$. The lower lines correspond to the expansion when one electron is removed, while the upper lines depict the contraction when an energy of  $10~\ensuremath{{\rm MeV}}$ is extracted. The lines correspond to post-bounce times of 73 ms (solid) and 240 ms (dashed). The labels "dump'' and "drain'' denote the density windows between the vertical, dotted lines, which roughly indicate where PNS convection deposits or extracts, respectively, lepton number and energy.

In Fig. C.1 we present the change of $\delta M$ when one electron or an energy of $10~\ensuremath{{\rm MeV}}$ is removed from the PNS at different radii (or densities, corresponding to different enclosed masses $M_{\rm p}$) for two representative times in the evolution of Model s15_32. Note that adding one electron or $10~\ensuremath{{\rm MeV}}$ would have the opposite effect. We see that removing this amount of energy always leads to contraction, while the extraction of one electron always causes expansion of the PNS. The latter finding can be understood from the fact that the electron is taken away but its degeneracy energy is assumed to remain in the PNS as thermal energy.

When discussing the effects of the convective transport of energy and lepton number on the PNS structure, it is interesting to compare the values of $\delta M$ in the drain region, where  $e_{\rm int}$ and  ${Y_{\rm lep}}$ drop, and in the dump region, where  $e_{\rm int}$ and  ${Y_{\rm lep}}$increase due to convection. For example, at $t=73{\rm ~ms}$, the removal of  $10~~\ensuremath{{\rm MeV}}$ from the drain region will lead to an increase of the mass enclosed by some large radius R by $(0.7{-}1.3)\times 10^{-24}$ g (Fig. C.1), while adding these $10~\ensuremath{{\rm MeV}}$ in the dump region will reduce $\delta M(R)$ by $(0.9{-}2.1)\times 10^{-24}$ g (Fig. C.1). This means a net decrease of $\delta M(R)$, and thus an expansion of the PNS. This behaviour is also seen at later times and in a similar way for the transport of ${Y_{\rm lep}}$.

   
Appendix D: Mixing scheme for reproducing PNS convection

In order to test the consequences of PNS convection on the PNS and supernova evolution in 1D simulations, we have developed a simple numerical "mixing algorithm'', which reproduces the energy and lepton number transport found in 2D models to a high degree.

In this algorithm we assume that any convective activity in the PNS leads instantaneously to a redistribution of energy and lepton number in the unstable region so that the convectively unstable gradients disappear, i.e.  $C_{\rm QL} = 0$ is established. Setting $\beta_{\rm diff}=1$and assuming ${\rm d}\left<{Y_{\rm lep}}\right>~/{\rm d}r~ = {\rm d}{Y_{\rm lep}}~/{\rm d}r~$ in Eq. (7), we find that the instability is mainly driven by a negative gradient of the entropy (including neutrino entropy). $C_{\rm QL} = 0$ then corresponds to a flat entropy profile, consistent with what is observed in the PNS convective region in 2D simulations, see Figs. 9 and 10. The numerical mixing scheme redistributes the energy in an unstable layer in a conservative way such that the entropy profile develops a plateau.

As a second constraint we introduce the empirically derived relation

$${Y_{\rm lep}}^{\rm new} = N {Y_{\rm lep}}^{\rm old} \left[1+ ... ... new} - e_{\rm int}^{\rm old}}{e_{\rm int}^{\rm old}}\right] ,$$ (D.1)

where "old'' denotes the lepton number ${Y_{\rm lep}}$ and internal energy  $e_{\rm int}$ before the application of the mixing scheme in a time step of the simulation, "new'' denotes the variables after such an application, and N is a global normalization factor which ensures that lepton number is globally conserved.

The mixing algorithm is executed in an operator splitting step after each hydrodynamic step. The algorithm detects regions with a negative entropy ($s+s_\nu$) gradient inside the PNS. In these regions energy and lepton number are then redistributed in such a way that Eq. (D.1), ${\rm d}(s+s_\nu~)/{\rm d}r~ \equiv 0$, and global energy and lepton number conservation are fulfilled (GR gravitational effects are ignored).

We find empirically that for $\alpha=2$ the evolution of the profiles of ${Y_{\rm lep}}$, $s+s_\nu$, and $e_{\rm int}$ inside the PNS for simulations with the mixing scheme reproduce very well the results of the 2D model, see Fig. 10. Near the boundaries of the convectively unstable layer small differences can appear. The over- and undershooting of dynamically moving matter into convectively stable layers can, of course, not be properly accounted for by this simple scheme. Another deficiency of the scheme is best visible shortly after the onset of convection around 50 ms after bounce, when the Brunt-Vaisälä frequency can reach values of up to 2 ms^-1 (Figs. 9 and D.1, panel b), and the convective layer has not yet fully developed and the transient starting phase with growing perturbations is not yet over. In this nonstationary, early phase of PNS convection the results with the mixing scheme cannot reproduce the ${Y_{\rm lep}}$-profiles of the 2D simulations very well (Fig. D.1, panel a); lepton number is transported more efficiently in the 2D models.

Figure D.1: Panel a) depicts the ${Y_{\rm lep}}$ profiles for three times around the onset of PNS convection, comparing a 1D model (s15s7b2), a 2D model (s15_32, for which angular-averaged values are plotted), and a 1D model with PNS mixing scheme (s15_mix). Panel b) shows  ${Y_{\rm lep}}$, the total entropy $s+s_\nu$, the radial velocity vr, and the Brunt-Väisälä-frequency $\omega _{\rm BV}$ (clockwise from top left to bottom left) in the PNS convective region for Model s15_32 shortly after the onset of PNS convection.

   
Appendix E: Perturbation growth during collapse

Figure E.1: Maximum fluctuations of the density within the core of $1.5~\ensuremath{M_{\odot}}$ at four times during the collapse of Model s15_64_p ( top). The initially imposed random perturbations are damped strongly during the subsonic collapse in the inner part of the Fe core. The corresponding density and velocity profiles (averaged in latitude) are shown in the middle and bottom panels, respectively. The short solid lines cutting the profiles in the bottom plot indicate the position of the sonic point.

Figure E.2: Standard deviations of the density fluctuations in lateral direction, $\sigma _\rho$, as defined in Buras et al. (2006, Eq. (24)) for different mass shells versus the (time-dependent) radius rM of these mass shells. The dashed lines represent fits according to $\sigma _\rho \propto r_M^{-1/2}$, the thick black bars mark the phases of the collapse when the infall velocity of a mass shell is more than 60% of the local sound speed, and the white bars mark the phases when the local infall velocity is supersonic.

The majority of our models were computed through the phases of core collapse, bounce, shock formation, and shock stagnation in spherical symmetry, and then continued in 2D with random perturbations added during the mapping from the 1D to the 2D grid. In contrast, we followed the evolution of Models s15_64_p and s15_64_r in two dimensions also during core collapse, shock formation, and early shock propagation. In this case the spherically symmetric density structure of the progenitor core was perturbed in each cell of the computational grid in a random way, in case of Model s15_64_p with an amplitude of 2%, in Model s15_64_r with 1% amplitude (Table 1). We will discuss here briefly the evolution of the perturbations during the core collapse phase, concentrating on Model s15_64_p. Since the initial spin is too low, rotation hardly affects the pre-bounce evolution of Model s15_64_r, and the latter model does not show any important rotation-induced differences with respect to the perturbation growth. When running the simulations continuously in 2D instead of mapping 1D post-bounce models on a 2D grid some milliseconds after bounce, we detect a slightly earlier onset of convection, which is seeded by the density and velocity fluctuations in the convectively unstable layers that develop below the neutrinosphere and in the neutrino-heating layer behind the stalled shock. Despite the earlier start of convection behind the shock, no significant differences of the evolution could be discovered between Models s15_32 and s15_64_p (Figs. 11-17). We therefore have included the results of the latter model in the relevant plots in Sect. 3, but have not specifically discussed them in the text.