Impacts of Collective Neutrino Oscillations on Supernova Explosions
Yudai Suwa, Kei Kotake, Tomoya Takiwaki, Matthias Liebendoerfer, Katsuhiko Sato
aa r X i v : . [ a s t r o - ph . H E ] J un Draft version August 8, 2018
Preprint typeset using L A TEX style emulateapj v. 8/13/10
IMPACTS OF COLLECTIVE NEUTRINO OSCILLATIONS ON CORE-COLLAPSE SUPERNOVA EXPLOSIONS
Yudai Suwa , Kei Kotake , Tomoya Takiwaki , Matthias Liebend¨orfer , and Katsuhiko Sato Draft version August 8, 2018
ABSTRACTBy performing a series of one- and two-dimensional (1-, 2D) hydrodynamic simulations with spec-tral neutrino transport, we study possible impacts of collective neutrino oscillations on the dynamicsof core-collapse supernovae. To model the spectral swapping which is one of the possible outcomeof the collective neutrino oscillations, we parametrize the onset time when the spectral swap begins,the radius where the spectral swap occurs, and the threshold energy above which the spectral inter-change between heavy-lepton neutrinos and electron/anti-electron neutrinos takes place, respectively.By doing so, we systematically study how the neutrino heating enhanced by the spectral swappingcould affect the shock evolution as well as the matter ejection. We also investigate the progenitordependence using a suite of progenitor models (13, 15, 20, and 25 M ⊙ ). We find that there is acritical heating rate induced by the spectral swapping to trigger explosions, which significantly differsbetween the progenitors. The critical heating rate is generally smaller for 2D than 1D due to themultidimensionality that enhances the neutrino heating efficiency. For the progenitors employed inthis paper, the final remnant masses are estimated to range in 1.1-1.5 M ⊙ . For our 2D model of the15 M ⊙ progenitor, we find a set of the oscillation parameters that could account for strong supernovaexplosions ( ∼ erg), simultaneously leaving behind the remnant mass close to ∼ . M ⊙ . Subject headings: hydrodynamics — neutrinos — radiative transfer — supernovae: general INTRODUCTION
Although the explosion mechanism of core-collapsesupernovae is not completely understood yet, currentmulti-dimensional (multi-D) simulations based on re-fined numerical models show several promising scenar-ios. Among the candidates are the neutrino heat-ing mechanism aided by convection and standing ac-cretion shock instability (SASI) (e.g., Marek & Janka2009; Bruenn et al. 2009; Suwa et al. 2010), the acous-tic mechanism (Burrows et al. 2007b), or the magnetohy-drodynamic (MHD) mechanism (e.g., Kotake et al. 2004,2006; Obergaulinger et al. 2006; Burrows et al. 2007a;Takiwaki et al. 2009). Probably the best-studied one isthe neutrino heating mechanism, whose basic conceptwas first proposed by Colgate & White (1966), and laterreinforced by Bethe & Wilson (1985) to take a currentlyprevailing delayed form.An important lesson from the multi-D simula-tions mentioned above is that hydrodynamic mo-tions associated with convective overturn (Herant et al.1994; Burrows et al. 1995; Janka & Mueller 1996;Fryer & Warren 2002, 2004) as well as the SASI (e.g.,Blondin et al. 2003; Scheck et al. 2006; Ohnishi et al.2006; Foglizzo et al. 2007; Murphy & Burrows 2008;Iwakami et al. 2008; Guilet et al. 2010; Fern´andez [email protected] Yukawa Institute for Theoretical Physics, Kyoto University,Oiwake-cho, Kitashirakawa, Sakyo-ku, Kyoto, 606-8502, Japan Division of Theoretical Astronomy, National AstronomicalObservatory of Japan, Mitaka, Tokyo 181-8588, Japan Center for Computational Astrophysics, National Astronom-ical Observatory of Japan, Mitaka, Tokyo 181-8588, Japan Department of Physics, University of Basel, Klingelbergstr.82, CH-4056 Basel, Switzerland The Institute for the Physics and Mathematics of the Uni-verse, the University of Tokyo, Kashiwa, Chiba, 277-8568, Japan National Institutes of Natural Sciences, Kamiyacho CentralPlace 2F, 4-3-13 Toranomon, Minato-ku, Tokyo, 105-0001, Japan . M ⊙ ) progen-itor of Woosley et al. (2002), and then for a 15 M ⊙ progenitor of Woosley & Weaver (1995) with a moder-ately rapid rotation imposed (Marek & Janka 2009). Byimplementing a multi-group flux-limited diffusion algo-rithm to the CHIMERA code (e.g., Bruenn et al. 2009),Yakunin et al. (2010) obtained explosions for a non-rotating 12 M ⊙ and 25 M ⊙ progenitor of Woosley et al.(2002). More recently, Suwa et al. (2010) pointed outthat a stronger explosion is obtained for a rapidlyrotating 13 M ⊙ progenitor of Nomoto & Hashimoto(1988) compared to the corresponding non-rotatingmodel, in which the isotropic diffusion source approx-imation (IDSA) for the spectral neutrino transport(Liebend¨orfer et al. 2009) is implemented in the ZEUScode.However, this success opens further new questions.First of all, the explosion energies obtained in these sim-ulations are typically underpowered by one or two ordersof magnitudes to explain the canonical supernova kineticenergy ( ∼ erg). Moreover, the softer nuclear equa-tion of state (EOS), such as of the Lattimer & Swesty(1991) (LS) EOS with an incompressibility K = 180 MeV SUWA ET AL.at nuclear densities is employed in those simulations. Ontop of evidence that favors a stiffer EOS based on nu-clear experimental data (Shlomo et al. 2006), the softEOS may not account for the recently observed massiveneutron star of ∼ M ⊙ (Demorest et al. 2010) (see themaximum mass for the LS180 EOS in O’Connor & Ott2011). With a stiffer EOS, the explosion energy may beeven lower as inferred from Marek & Janka (2009) whodid not obtain the neutrino-driven explosion for theirmodel with K = 263 MeV. What is then missing fur-thermore? We may get the answer by going to 3D simu-lations (Nordhaus et al. 2010) or by taking into accountnew ingredients, such as exotic physics in the core of theprotoneutron star (Sagert et al. 2009), viscous heating bythe magnetorotational instability (Thompson et al. 2005;Masada et al. 2011), or energy dissipation via Alfv´enwaves (Suzuki et al. 2008).Joining in these efforts, we explore in this study thepossible impacts of collective neutrino oscillations on en-ergizing the neutrino-driven explosions. The collectiveneutrino oscillations, i.e. neutrinos of all energies thatoscillate almost in phase, are attracting great attention,because they can induce dramatic observable effects suchas a spectral split or swap (e.g., Raffelt & Smirnov 2007;Duan et al. 2008; Dasgupta et al. 2008, and see refer-ences therein). They are predicted to emerge as a distinctfeature in their energy spectra (see Duan et al. 2010;Dasgupta 2010, for reviews of the rapidly growing re-search field and collective references therein). Amonga number of important effects possibly created by theself-interaction, we choose to consider the effect of spec-tral splits between electron- ( ν e ), anti-electron neutrinos(¯ ν e ), and heavy lepton neutrinos ( ν x , i.e., ν µ , ν τ and theiranti-particles) above a threshold energy (e.g., Fogli et al.(2007)). Since ν x ’s have higher average energies than theother species in the postbounce phase, the neutrino fla-vor mixing would increase the effective energies of ν e and¯ ν e , and hence increase the neutrino heating rates in thegain region. A formalism to treat the neutrino oscillationin the Boltzmann neutrino transport is given in Yamada(2000); Strack & Burrows (2005), but difficult to imple-ment. To just mimic the effects in this study, we performthe spectral swap by hand as a first step. By changingthe average neutrino energy, h ǫ ν x i , as well as the positionof the neutrino spheres ( R ν x ) in a parametric manner,we hope to constrain the parameter regions spanned by h ǫ ν x i and R ν x in which the additional heating given bythe collective neutrino oscillations could have impacts onthe explosion dynamics. Our strategy is as follows. Byperforming a number of 1D simulations, we will firstlyconstrain the parameter regions to some extent. Here wealso investigate the progenitor dependence using a suiteof progenitor models (13, 15, 20, and 25 M ⊙ ). Aftersqueezing the condition in the 1D computations, we in-clude the flavor conversions in 2D simulations to see theirimpacts on the dynamics, and we also discuss how thecritical condition for the collective effects in 1D can besubject to change in 2D.The paper opens with descriptions of the initial mod-els and the numerical methods focusing how to model thecollective neutrino oscillations (Section 2). The main re-sults are shown in Section 3. We summarize our resultsand discuss their implications in Section 4. NUMERICAL METHODS
Hydrodynamics
The employed numerical methods are essentially thesame as those in our previous paper (Suwa et al. 2010).For later convenience, we briefly summarize them in thefollowing. The basic evolution equations are written as,d ρ d t + ρ ∇ · v = 0 , (1) ρ d v d t = −∇ P − ρ ∇ Φ , (2)d e ∗ d t + ∇ · [( e ∗ + P ) v ] = − ρ v · ∇ Φ + Q E , (3)d Y e d t = Q N , (4) △ Φ = 4 πGρ, (5)where ρ, v , P, v , e ∗ , Φ, are density, fluid velocity, gas pres-sure including the radiation pressure of neutrinos, totalenergy density, gravitational potential, respectively. Thetime derivatives are Lagrangian. As for the hydro solver,we employ the ZEUS-2D code (Stone & Norman 1992)which has been modified for core-collapse simulations(e.g., Suwa et al. 2007b,a, 2009; Takiwaki et al. 2009). Q E and Q N (in Equations (3) and (4)) represent thechange of energy and electron fraction ( Y e ) due to theinteractions with neutrinos. To estimate these quanti-ties, we implement spectral neutrino transport using theisotropic diffusion source approximation (IDSA) scheme(Liebend¨orfer et al. 2009). The IDSA scheme splits theneutrino distribution into two components, both of whichare solved using separate numerical techniques. We ap-ply the so-called ray-by-ray approach in which the neu-trino transport is solved along a given radial direction as-suming that the hydrodynamic medium for the directionis spherically symmetric. Although the current IDSAscheme does not yet include ν x and the inelastic neu-trino scattering with electrons, these simplifications savea significant amount of computational time compared tothe canonical Boltzmann solvers (see Liebend¨orfer et al.(2009) for more details). Following the prescription inM¨uller et al. (2010), we improve the accuracy of the to-tal energy conservation by using a conservation form inequation (3), instead of solving the evolution of internalenergy as originally designed in the ZEUS code. Numer-ical tests are presented in Appendix A.The simulations are performed on a grid of 300 log-arithmically spaced radial zones from the center up to5000 km and 128 equidistant angular zones covering0 ≤ θ ≤ π for two-dimensional simulations. For the spec-tral transport, we use 20 logarithmically spaced energybins reaching from 3 to 300 MeV. Spectral swapping
As mentioned in §
1, we introduce a spectral inter-change from heavy-lepton neutrinos ( ν µ , ν τ and theirantineutrinos, collectively referred as ν x hereafter) toelectron-type neutrinos and antineutrinos, namely ν x → ν e and ¯ ν x → ¯ ν e . Instead of solving the transportmpacts of collective neutrino oscillations on supernova explosions 3equations for ν x , we employ the so-called light-bulb ap-proximation and focus on the optically thin region out-side the neutrinosophere (e.g., Janka & Mueller 1996;Ohnishi et al. 2006).According to Duan et al. (2010), we set the thresholdenergy, ǫ th , to be 9 MeV, above which the spectral swaptakes place. Below the threshold, the neutrino heating isestimated by the spectral transport via the IDSA scheme.Above the threshold, the heating rate is replaced by Q E ∝ Z ∞ ǫ th dǫ ν ǫ [ j ( ǫ ν ) + χ ( ǫ ν )] f ν ( r, ǫ ν ) , (6)where j and χ are the neutrino emissivity and absorp-tivity, respectively, and f ν ( r, ǫ ν ) corresponds to the neu-trino distribution function for ν x with ǫ ν being energiesof electron neutrinos and antineutrinos. In the light-bulbapproach, it is often approximated by the Fermi-Diracdistribution with a vanishing chemical potential (e.g.,Ohnishi et al. 2006) as, f ν ( r, ǫ ν ) = 1 e ǫ ν /kT νx + 1 g ( r ) , (7)where k , T ν x are the Boltzmann constatn and the neu-trino temperature, respectively. g ( r ) is the geometricfactor, g ( r ) = 1 − (cid:2) − ( R ν x /r ) (cid:3) / which is taken intoaccount for the normalization, with R ν x being the radiusof the neutrinosphere. The neutrino luminosity of ν x atthe infinity is the given as L ν x = 2 . × (cid:18) h ǫ ν x i
15 MeV (cid:19) (cid:18) R ν x
30 km (cid:19) erg s − , (8)where h ǫ ν x i = R ∞ dǫ ν x ǫ ν x f ν ( ǫ ν x ) / R ∞ dǫ ν x ǫ ν x f ν ( ǫ ν x ) isthe average energy of emitted neutrinos. The positionwhere the spectral swapping sets in is fixed at 100 km(around the gain radius) and the onset time is varied asa parameter, t s =100, 200, and 300 ms after bounce.In fact, the threshold energy depends on the neu-trino luminosities, spectra and oscillation parameters(see, e.g., Duan et al. 2010, and references therein) withconserved net ν e flux (i.e., the lepton number conserva-tion). However, the conservation of lepton number is toocomplicated to satisfy in the dynamical simulation be-cause the neutrino spectrum and the luminosity evolvewith time. In order to focus on the hydrodynamic fea-tures affected by the spectral modulation induced by theswapping, we simplify just a single threshold energy inthis work.To summarize, the parameters that we use to mimicthe spectral swapping are the following three items, (i) R ν x which is the radius of the neutrinosphere of ν x , (ii) h ǫ ν x i which is the average energy of ν x , and (iii) t s whichis the time when the spectral swapping sets in. RESULT
One-dimensional models
1D without spectral swapping
In this subsection, we first outline the 1D collapse dy-namics without spectral swapping. We take a 13 M ⊙ progenitor (Nomoto & Hashimoto 1988) as a reference.At around 112 ms after the onset of gravitational col-lapse, the bounce shock forms at a radius of ∼
10 km < ε ν > [ M e V ] Time after Bounce [ms] 0 0.2 0.4 0.6 0.8 1 L ν [ e r g s - ] ν e anti- ν e Figure 1.
Time evolution of the neutrino luminosity (top panel)and average energy (bottom panel) for ν e (red-solid line) and ¯ ν e (blue-dashed line). with an enclosed mass of ∼ . M ⊙ . The central den-sity at this time is ρ c = 3 . × g cm − . The shockpropagates outwards but finally stalls at a radius of ∼
100 km. Due to the decreasing accretion rate throughthe stalled shock, the shock can be still pushed outward.However, after some time, the shock radius begins toshrink. The ratio of the advection timescale, τ adv , andthe heating timescale, τ heat , is an important indicatorfor the criteria of neutrino driven explosion (Buras et al.2006; Marek & Janka 2009; Suwa et al. 2010). In our 1Dsimulations, τ adv /τ heat is generally smaller than unity inthe postbounce phase. This is the reason why our 1Dsimulations do not yield a delayed explosion. This alsothe case for the other progenitors (15, 20, and 25 M ⊙ ) in-vestigated in this study. As for the accretion phase (laterthan ∼
50 ms after the bounce), the typical neutrino lu-minosity at r = 5000 km is 3 × erg s − for both ν e and ¯ ν e , and the typical average energy is h ǫ ν e i ≈ h ǫ ¯ ν e i ≈
12 MeV as shown in Figure 1. Figure 2indicates the resultant neutrino luminosity spectrum at100 ms after the bounce.
1D with spectral swapping
The investigated models with the spectral swappingare summarized in Table 1. As already mentioned, themodel parameters are the neutrinosphere radius ( R ν x ),the average energy of neutrinos ( h ǫ ν x i ), and the onsettime of the spectral swapping ( t s ). The model names in-clude these parameters; “NH13” represents the progen-itor model, “R..” represents R ν x in units of km, “E..”represents h ǫ ν x i in MeV, “T..” represents t s in ms, andthe last letter “S” represents 1D (spherical symmetry).Figure 3 presents the time evolution of the mass shellsfor models NH13R30E12T100S and NH13R30E13T100S.The difference between these panels is the average en-ergies of neutrinos, h ǫ ν x i = 12 MeV for the top paneland 13 MeV for the bottom panel. The thick solid lines Note that 0 . M ⊙ is rather high value that is due to ap-proximations employed in our simulation. We omit the electronscattering by neutrinos and general-relativistic effects, which leadsmaller inner core mass at bounce (see Liebend¨orfer et al. 2001;Thompson et al. 2003). In addition, more improved electron cap-ture treatment would lead even smaller (Langanke et al. 2003). SUWA ET AL. d L ν / d ε ν [ e r g M e V - ] ε ν [MeV]< ε ν >=15MeV ν e anti- ν e ν x with Figure 2.
The neutrino luminosity spectrum of ν e (red-solid line)and ¯ ν e (blue-dashed line) without the spectral swapping at 100 msafter the bounce. For comparison, we show the injected luminosityspectrum of ν x with h ǫ ν x i = 15MeV, which will be swapped withthe original spectrum of ν e and ¯ ν e at ǫ ν > represent the radial position of shock waves. Regardlessof a small difference of h ǫ ν x i , model NH13R30E13T100Sshows a shock expansion after the manual spectral swap-ping is switched on (see the thick line in the bottompanel of Figure 3), while the stalled shock does not re-vive for model NH13R30E12T100S (top panel). Thissuggests that there is a critical condition for the suc-cessful explosion induced by the spectral swapping. Inthe bottom panel, the regions enclosing the mass of M r ∼ . M ⊙ (thin black line) corresponds to the so-called mass cut, which could be interpreted as the finalmass of the remnant. The fact that a clear mass cutemerges in model NH13R30E13T100S indicates that aneutron star will be left behind in this model. Such a def-inite mass-cut has been observed in Kitaura et al. (2006)who reported a successful neutrino-driven explosion (in1D) for a lighter progenitor star, which is, however, diffi-cult to realize for more massive stars in 2D (e.g., Figure2 in Marek & Janka (2009) and Figure 1 in Suwa et al.(2010)).As a tool to measure the strength of an explosion, wedefine a diagnostic energy that refers to E diag = Z D dV (cid:18) ρ | v | + e − ρ Φ (cid:19) , (9)where e is internal energy, D represents the domain inwhich the integrand is positive. Figure 4 shows the timeevolution of E diag for some selected models. The diag-nostic energy increases with time for the green-dottedline, which turns to decrease for the red line, noting thatthe difference between the pair of models is ∆ h ǫ ν x i = 1MeV. The blue-dashed line (model NH13R30E15T100S)has h ǫ ν x i = 15 MeV and reaches larger E diag than thegreen line (NH13R30E13T100S; h ǫ ν x i = 13 MeV). On theother hand, the later injection of the spectral swappingleads to smaller E diag , i.e. the brown-dot-dashed line( t s =200 ms) shows smaller E diag than the blue-dashedline ( t s =100 ms). For models that experience earlierspectral swapping with higher neutrino energy, the diag-nostic energy becomes higher in an earlier stage, as it isexpected.
10 100 1000 10000 100 150 200 250 300 350 400 R a d i u s [ k m ] Time after boune [ms]N13R30E12T100 10 100 1000 10000 100 150 200 250 300 350 400 R a d i u s [ k m ] Time after boune [ms]N13R30E13T100
Figure 3.
Time evolution of mass shells for NH13R30E12T100S(top) and NH13R30E13T100S (bottom). The Black thin line cor-responds to 1 . M ⊙ and the black thick line represents the shockwave position, respectively. The difference between these panelsis the average energies of neutrinos, h ǫ ν x i = 12 MeV for the toppanel and 13 MeV for the bottom panel. Looking at Figure 4 again, E diag for the explodingmodels seems to show a saturation with time. Thesecurves can be fitted by the following function, E diag ( t ) = E ∞ diag (1 − e − at + b ) , (10)where E ∞ diag is a converging value of E diag , a and b are the fitting parameters. As for NH13R30E13T100S, E ∞ diag = 8 . × erg. This fitting formula allows usto estimate the final diagnostic energy especially for thestrongly exploding models whose diagnostic energy wecannot estimate in principle because the shock goes be-yond the computational domains ( r < R ν x and h ǫ ν x i (equation (8)). The gray lines correspond to the neutrinoluminosities determined by the pairs of R ν x and h ǫ ν x i which is 1 to 5 × erg s − from bottom to top lines.Circles and crosses correspond to the exploding and non-exploding models, respectively. Not surprisingly, explo-sions are more easier to be obtained for higher neutrinoluminosity.As is well known, the combination of h ǫ ν x i and L ν x ismpacts of collective neutrino oscillations on supernova explosions 5 -4 -3 -2
100 150 200 250 300 350Time after bounce [ms]0.511.522.53 D i a gno s ti c E n e r gy [ e r g ] N13R30E12T100SN13R30E13T100SN13R30E15T100SN13R30E15T200S
Figure 4.
Diagnostic energies as functions of time. Red-solid, green-dotted, blue-dashed and brown-dot-dashed linescorrespond to models NH13R30E12T100S, NH13R30E13T100S,NH13R30E15T100S, and NH13R30E15T200S, respectively. The ν x average energies of three lines other than brown-dot-dashed linediffer from each other, i.e. h ǫ ν x i = 12 MeV (red solid), 13 MeV(green dotted), and 15 MeV (blue dashed), respectively. As forthe brown-dot-dashed line, only the onset time of spectral swap-ping is different from the blue-dashed line. The red-solid line showsonly an oscillation, while the other lines show increasing diagnosticenergy.
10 15 20 25 0 10 20 30 40 50 60 < ε ν > [ M e V ] R ν [km]L ν =10 erg s -1 L ν =5 × erg s -1 ExplodingNonexploding
Figure 5.
Summary of 1D models with t s = 100 ms after thebounce. Circles indicate the exploding models, while crosses shownon-exploding models. Gray solid lines correspond to the luminos-ity of ν x calculated by Eq. (8), which are 1, 2, 3, 4, and 5 × erg s − from bottom to top. an important quantity to diagnose the success or failureof explosions, because the neutrino heating rate in theso-called gain region, Q + ν , is proportional to (cid:10) ǫ ν x (cid:11) L ν x (e.g., equation (23) in Janka (2001)).Figure 6 shows E ∞ diag as a function of h ǫ ν x i L ν x . Notein the plot that we set the horizontal axis not as (cid:10) ǫ ν x (cid:11) L ν x but as h ǫ ν x i L ν x so that we can deduce the following de-pendence more clearly and easily . In this figure, let usfirst focus on red pluses, green crosses, and blue squares (cid:10) ǫ ν x (cid:11) (cid:0) ≡ R ∞ dǫ ν x ǫ ν x f ν ( ǫ ν x ) / R ∞ dǫ ν x ǫ ν x f ν ( ǫ ν x ) (cid:1) and h ǫ ν x i can be simply connected as (cid:10) ǫ ν x (cid:11) = 2 . h ǫ ν x i for the neutrinospectrum of equation (7). E ∞ d i a g [ e r g ] < ε ν > L ν [10 MeV erg s -1 ]t s =100mst s =150mst s =200ms2D 0 0.5 1 1.5 2 2.5 0 1 2 3 4 5 6 7 E ∞ d i a g [ e r g ] < ε ν > L ν [10 MeV erg s -1 ] Figure 6.
Diagnostic energies for exploding models at severalhundred seconds after the bounce. Red points, green crosses, andblue squares correspond to models with t s =100, 150, and 200 ms,respectively. Red circles represent the result of 2D simulation (seetext for details). whose difference is characterized by t s (2D results (filledcircles) will be mentioned in the later section). Red( t s =100 ms), green ( t s =150 ms), and blue ( t s =200ms) points have a clear correlation with h ǫ ν x i L ν x . Or-ange and light-blue regions represent the non-explodingregions for red and blue points, respectively. Both ofthem show that the minimum E ∞ diag decreases with t s ,indicating that the critical values of h ǫ ν i L ν for explo-sion sharply depends on t s . This is because the massoutside the shock wave gets smaller with time so thatthe minimum energy to blow up star gets smaller too.By the same reason, E diag becomes larger as t s becomessmaller given the same h ǫ ν x i L ν x . To obtain a larger E ∞ diag , the earlier spectral swapping is more preferential.Figure 7 shows the neutrino heating rate and the den-sity distribution of NH13R30E13T100S for 10 ms and250 ms after t s (=100 ms after the bounce). As theshock wave propagates outward, the density in the gainregion sharply drops (e.g., 100-200km, dashed blue line),leading to the suppression of the heating rate (dashedred line). This is the reason of the saturation in E diag asshown in Figure 4.The remnant mass is an important indicator to diag-nose the consequences of the explosion in producing ei-ther a neutron star or a black hole. The last two linesin Table 1 show the integrated masses in the regions of ρ ≥ g cm − at t = t s and t = ∞ . The latter one isestimated by the fitting as M ( t ) = M ∞ (1 + e − ct + d ) , (11)where c and d are the fitting parameters. For theexploding models, M ∞ becomes generally smallerthan M t = t s because of the mass ejection. Exceptionsare weakly exploding models (NH13R20E15T150S,NH13R20E15T200S, NH13R30E13T100S, andNH13R50E11T100S), in which the mass accretioncontinues after t s and stops eventually at late time(maximum masses are presented in Table 1). Forthe nonexploding models, the remnant mass simplyincreases with time. Regarding the 13 M ⊙ progenitor SUWA ET AL. -2-1.5-1-0.5 0 0.5 1 1.5 2 100 1000 10 Q ν R [ e r g s - ] ρ [ g c m - ] Radius [km]Q ν R for 10 ms after t s Q ν R for 250 ms after t s ρ for 10 ms after t s ρ for 250 ms after t s Figure 7.
Heating rate at 10 ms (red-solid line) and 250 ms (red-dashed line) after the bounce for the model NH13R30E13T100Sand the density profile (blue-solid and dashed lines). As the densitydecreases due to the neutrino driven wind, the heating rate alsodecreases. ivestigated in this section, the remnant masses in modelsthat produce strong explosion ( E ∞ diag & erg), areconsiderably smaller (1.1-1.2 M ⊙ ) if compared to thetypical mass as of observed neutron stars ∼ . M ⊙ (Lattimer & Prakash 2007). This may simply reflect thelight iron core ( ∼ . M ⊙ ) inherent to the progenitoror the existence of mass accretion induced by thematter fallback after the explosion. Now we move onto investigate the progenitor dependence in the nextsection. The progenitor dependence
In addition to the 13 M ⊙ progenitor byNomoto & Hashimoto (1988), we are going to investigatethe progenitor dependence in 1D simulations. The com-puted models are NH15 (15 M ⊙ ) (Nomoto & Hashimoto1988), s15s7b2 (15 M ⊙ ) (Woosley & Weaver 1995),s15.0 (15 M ⊙ ), s20.0 (20 M ⊙ ), and s25.0 (25 M ⊙ )(Woosley et al. 2002), which are listed in Table 2. Thefirst sets of characters for these models indicate theprogenitors as, • NH: (Nomoto & Hashimoto 1988) • WW: (Woosley & Weaver 1995) • WHW: (Woosley et al. 2002)Figure 8 depicts density profiles of these progenitors100 ms after the bounce as a function of the enclosedmass ( M r ). It can be seen that the density profiles for M r . . M ⊙ are almost insensitive to the progenitormasses despite the difference in the pre-collapse phase(see, e.g., Figure 1 of Burrows et al. 2007b). On the otherhand, the profiles of M r & . M ⊙ differ between pro-genitors so that the critical heating rates and E ∞ diag areexpected to be different also. In Figure 8, the envelopeof WHW25 is shown to be thickest, while the envelopeof NH13 is thinnest.Figure 9 shows the critical heating rates as a functionof the progenitor masses. In agreement with intuition,the critical heating rate for models WHW25 and NH13belongs to the high and low ends, respectively. However,
0 0.5 1 1.5 2 2.5 D e n s it y [ g c m - ] Enclosed Mass [M ⊙ ] N13 (NH88)N15 (NH88)S15 (WW95)H15 (WHW02)H20 (WHW02)H25 (WHW02) Figure 8.
Density profiles of investigated progenitors 100 ms afterthe bounce as functions of the enclosed mass. < ε ν > L ν [ e r g s - M e V ] M [M ⊙ ] NH88WW95WHW02 Figure 9.
The critical heating rate, h ǫ ν x i L ν x , as a func-tion of the progenitor mass, M . Circles are progeni-tors from Nomoto & Hashimoto (1988), the square is fromWoosley & Weaver (1995), and crosses are from Woosley et al.(2002), respectively. The error bars represent the distance betweenthe last failing and the first exploding model in our grid of models.The symbols locate at the centers of error bars. The error bar issmall for model NH13 because we calculated a more refined gridof models for 13 M ⊙ progenitor (Table 1) than for the 15-25 M ⊙ progenitors (Table 2). the critical heating rate for model WHW20 is almost thesame as the one for model NH13 although the envelopeof model WHW20 is much thicker than model NH13 (seeFigure 8). Our results show that the critical heating rateis indeed affected by the envelope mass, however, therelation is not one-to-one. It is also interesting to notethat the critical heating rates for 15 M ⊙ progenitors ofWW15, WHW15 and NH15, are different by a factor of ∼
3, which may send us a clear message that the accurateknowledge of supernova progenitors is also pivotal to pindown the supernova mechanism.The integrated masses with ρ ≥ g cm − for t = t s and t = ∞ are listed in the last two lines in Table 2 andFigure 10. The tendencies are the same as found withNH13. As for model WHW25, we obtain results with E ∞ diag > erg and M ∞ > . M ⊙ , simultaneously. Two-dimensional models mpacts of collective neutrino oscillations on supernova explosions 7
Table 1
1D simulationsModel Dimension R ν h ǫ ν x i L ν t s Explosion E ∞ diag M t = t s M ∞ [km] [MeV] [10 erg s − ] [ms] [10 erg] [ M ⊙ ] [ M ⊙ ]NH13R10E15T100S 1D 10 15MeV 0.29 100 No — 1.18 —NH13R10E17T100S 1D 10 17MeV 0.48 100 No — 1.18 —NH13R10E18T100S 1D 10 18MeV 0.60 100 No — 1.18 —NH13R10E19T100S 1D 10 19MeV 0.75 100 Yes
Yes
Yes
Yes < . Yes < . Yes
Yes
Yes < . Yes
Yes
Yes
Yes
Yes
Yes
Yes
Yes < . Yes M N S [ M ⊙ ] M [M ⊙ ] NH88WW95WHW02 Figure 10.
The final NS masses as a function of the progeni-tor mass, M . Circles are progenitors from Nomoto & Hashimoto(1988), squares are from Woosley & Weaver (1995), and crosses arefrom Woosley et al. (2002), respectively. Here we discuss the effects of spectral swapping in 2D(axisymmetric) simulations. Since our 2D simulations,albeit utilizing the IDSA scheme, are still computation-ally expensive, it is not practicable to perform a system-atic survey in 2D as we have done in 1D simulations.Looking at Figure 9 again, we choose models WHW15 (Woosley et al. 2002) and NH13 (Nomoto & Hashimoto1988), whose critical heating rate belong to the high andlow ends, respectively.
2D without spectral swapping
The basic hydrodynamic picture is the same with 1Dbefore the shock-stall (e.g., till .
10 ms after bounce).After that, convection as well as SASI sets in betweenthe stalled shock and the gain radius, which leads tothe neutrino-heated shock revival for model NH13 (e.g.,Suwa et al. (2010)). While for model WHW15, the po-sition of the stalled shock, following several oscillations,begins to shrink at &
400 ms after bounce.Even after the shock revival, it should be emphasizedthat the shock propagation for model NH13 is the so-called “passive” one (Buras et al. 2006). This means thatthe amount of the mass ejection is smaller than the accre-tion in the post-shock region of the expanding shock (seemotions of mass shells in the post-shock region of Figure1 in Suwa et al. (2010)). Some regions have a positivelocal energy (Eq. (9)), but the volume integrated valueis quite as small as . erg at the maximum. In orderto reverse the passive shock into an active one it is mostimportant to energize the explosion in some way. Usingthese two progenitors that produce a very weak explo-sion (model NH13) and do not show even a shock revival(model WHW15), we hope to explore how the dynamics SUWA ET AL. Table 2
Progenitor dependenceModel Dimension R ν h ǫ ν x i L ν t s Explosion E ∞ diag M t = t s M ∞ [km] [MeV] [10 erg s − ] [ms] [10 erg] [ M ⊙ ] [ M ⊙ ]NH15R30E11T100S 1D 30 11MeV 0.76 100 No — 1.34 —NH15R30E12T100S 1D 30 12MeV 1.07 100 No — 1.34 —NH15R30E13T100S 1D 30 13MeV 1.48 100 Yes < . Yes
Yes
Yes
Yes
Yes
Yes
Yes
Yes
Yes < . Yes
Yes would change when the spectral swapping is switched on.
2D with spectral swapping
Table 3 shows a summary for our 2D models, in whichthe last character of each model (A) indicates “Axisym-metry”. Models NH13A and WHW15A are 2D modelswithout spectral swapping for NH13 and WHW15, re-spectively.As in 1D, the onset of the spectral swapping is takento be t s = 100 ms after bounce. At this time, modelNH13 shows the onset of the gradual shock expansionwith a small diagnostic energy of E diag ∼ × erg,and the shock radius is located at ∼
300 km. As formodel WHW15, there is no region with a positive localenergy (e.g., Eq. (9)) and the shock radius is ∼
200 km.The density profile for this model is essentially same asthe one in the 1D counterpart (see Figure 8) but withsmall angular density modulations due to convection.In Figure 6, red filled circles represent E ∞ diag for modelNH13. It can be seen that the critical heating rate toobtain E ∞ diag ∼ erg is smaller for 2D than the corre-sponding 1D counterparts (compare the heating rates for h ǫ ν x i L ν x ∼ . × MeV erg s − ). In fact, modelswith h ǫ ν x i L ν x . . × MeV erg s − fail to explodein 1D, but succeed in 2D (albeit with a relatively small E ∞ diag less than 10 erg). As opposed to 1D, it is ratherdifficult in 2D to determine a critical heating rate due tothe stochastic nature of the explosion triggered by SASIand convection. In our limited set of 2D models, the crit-ical heating rate is expected to be close to h ǫ ν x i L ν x ∼ . × MeV erg s − , below which the shock does notrevive (e.g., h ǫ ν x i L ν x . × MeV erg s − , is thelowest end in the horizontal axis in the figure).As seen from Figure 6, E ∞ diag becomes visibly larger for 2D than 1D especially for a smaller h ǫ ν x i L ν x . Asthe heating rates become larger, the difference between1D and 2D becomes smaller because the shock revivaloccurs almost in a spherically symmetric way (beforeSASI and convection develop non-linearly). In Table 3,it is interesting to note that model NH13R30E11T100Afails to explode, while we observed the shock-revival forthe corresponding model without the spectral swapping(model NH13A). This is because the heating rate ofmodel NH13R30E11T100A is smaller than NH13A dueto the small h ǫ ν x i , which can make it more difficult totrigger ν x explosions. On the other hand, if the energygain due to the swap is high enough (i.e., for models withgreater than E12 in Table 3), the swap can facilitate ex-plosions.Figure 11 depicts the entropy distributions for mod-els NH13A (top panel) and NH13R30E13T100A (bot-tom panel). It can be seen that model NH13A showsa unipolar-like explosion (see also Suwa et al. 2010),while model NH13R30E13A explodes rather in a spher-ical manner as mentioned above. Model NH13A expe-riences several oscillations aided by SASI and convec-tion before explosion, while the stalled shock for modelNH13R30E13T100A, turns into expansion shortly afterthe onset of the spectral swapping. In fact, the shapesof hot bubbles behind the expanding shock are shown tobe barely changing with time (bottom panel), which in-dicates a quasi-homologous expansion of material behindthe revived shock.Figure 12 shows the time evolution of massshells for models NH13A (thin-gray lines) andNH13R30E13T100A (thin-orange lines). Black andred thick lines represent the shock position at the northpole for each models. The mass shells for model NH13Acontinue to accrete to the PNS, since the shock passively expands outwards as already mentioned. Due to thismpacts of collective neutrino oscillations on supernova explosions 9 Figure 11.
Time evolution of the entropy distributions.
Top : NH13A without the spectral swapping for 100, 200, 300, and 450 ms afterbounce from left to right.
Bottom : NH13R30E13T100A for 100, 150, 200, 250 ms after bounce (corresponding to 0, 50, 100, 150 ms afterthe onset of the spectral swapping.) continuing mass accretion, the remnant for this modelwould be a black hole instead of a neutron star. Onthe other hand, model NH13R30E13T100A shows amass ejection with a definite outgoing momentum inthe postshock region so that the remnant could bea neutron star. Unfortunately however, we cannotpredict the final outcome due to the limited simulationtime. A long-term simulation recently done in 1D (e.g.,Fischer et al. (2010)) should be indispensable also forour 2D case. This is, however, beyond the scope of thispaper. Here let us discuss a validity of the parameters for thespectral swap that we have assumed so far. For example,the criteria of explosion for model NH13R30E12T100Awas L ν x ≈ . × erg s − and h ǫ ν x i ≈
12 MeV. Thesevalues are even smaller than the typical values obtainedin 1D Boltzmann simulations (e.g., Liebend¨orfer et al.(2004)), which show L ν x ≈ × erg s − and q(cid:10) ǫ ν x (cid:11) ≈
20 MeV (i.e. h ǫ ν x i ≈
14 MeV with a vanishing chemicalpotential) earlier in the postbounce phase. Therefore thespectral swapping, if it would work as we have assumed,0 SUWA ET AL.
10 100 1000 0 50 100 150 200 250 300 350 400 R a d i u s [ k m ] Time after boune [ms]
Figure 12.
Time evolution of mass shells for NH13A (thin-graylines) and NH13R30E13T100A (thin-orange lines). Black and redthick lines represent the shock position at the north pole. may be a potential to assist explosions.It should be noted that the critical heating rate in thisstudy might be too small due to the approximation ofthe light-bulb scheme. In this scheme, we can include thegeometrical effect of the finite size of the neturinosphereas in Eq. (7), but can not include the back reactionby the matter, i.e. the absorption of neutrino. Somefraction of neutrinos, in fact, are absorbed in the gainregion and the neutrino luminosity decreases with theradius. We omit this effect in this study so that theheating rate might be overestimated in the simulationwith the spectral swapping. Thus, the fully consistentsimulation including the spectral swapping is necessaryfor more realistic critical heating rate, which is beyondthe scope of this study.Finally we discuss the 15 M ⊙ progenitor labeled byWHW15. As mentioned, this progenitor fails to explodewithout spectral swapping even in 2D . Figure 13 showsthe entropy distributions of WHW15A (left; nonexplod-ing) and WHW15R30E15T100A (right; exploding) for220 ms after the bounce (corresponding to 120 ms after t s for model WHW15R30E15T100A). The model with R ν x = 30 km and h ǫ ν x i = 14 MeV does not explodein 1D but explodes in 2D (compare Table 2 and 3).Again, the mulitidimensionality helps the onset of ex-plosion. The critical heating rate in 2D is in the rangeof 2 . ≤ h ǫ ν x i L ν x / (10 MeV erg s − ) ≤ .
9, whileit is 3 . ≤ h ǫ ν x i L ν x / (10 MeV erg s − ) ≤ . ∼ ν x luminosity and average energy to obtain explosion are L ν x ∼ × erg s − and h ǫ ν x i ∼
14 MeV (correspond-ing to q(cid:10) ǫ ν x (cid:11) ∼
20 MeV), which are close to the resultsobtained in a 1D Boltzmann simulation (Sumiyoshi et al.2005) for a 15 M ⊙ progenitor . The diagnostic energy This is consistent with a very recent result byObergaulinger & Janka (2011), who performed 2D simulations ofmodel WHW15 with spectral neutrino transport Note that the progenitor employed in Sumiyoshi et al. (2005)is WW95, so that the direct comparison may not be fair. However,the critical heating rate in 1D for WW15 is smaller than WHW15(Figure 9) and the mass of the envelope is thicker for WHW15 than
Figure 13.
The entropy distributions of WHW15A (left) andWHW15R30E15T100A (right) for 220 ms after the bounce. as well as the estimated remnant masses are listed in thelast three columns in Table 3. E ∞ diag (as well as M ∞ diag ) forexploding models is shown to be larger than the modelseries of NH13. As a result, some of the 2D models forWHW15 produce strong explosions ( E ∞ diag ∼ erg),while simultaneously leaving behind a remnant of 1.34–1.52 M ⊙ . We think that it is only a solution acciden-tally found by our parametric explosion models. Howeveragain, the critical heating rates that require to assist theneutrino-driven explosion via the spectral swapping arenever far away from the ones obtained in the Boltzmannsimulations. We hope that our exploratory results maygive a momentum to supernova theorists to elucidate theeffects of collective neutrino oscillations in a more con-sistent manner. SUMMARY AND DISCUSSION
We performed a series of one- and two-dimensionalhydrodynamic simulations of core-collapse supernovaewith spectral neutrino transport via the IDSA scheme.To model the spectral swapping which is one of thepossible outcomes of the collective neutrino oscilla-tions, we parametrized the onset time when the spec-tral swap begins, the radius where the spectral swaptakes place, and the threshold energy above which thespectral interchange between heavy-lepton neutrinos andelectron/anti-electron neutrinos occurs. By doing so, wesystematically studied the shock evolution and the mat-ter ejection due to the neutrino heating enhanced byspectral swapping. We also investigated the progenitordependence using a suite of progenitor models (13, 15,20, and 25 M ⊙ ). With these computations, we foundthat there is a critical heating rate induced by the spec-tral swapping to trigger explosions, which differs betweenthe progenitors. The critical heating rate is generallysmaller for 2D than 1D due to the multidimensionalitythat enhances the neutrino heating efficiency (see alsoJanka & Mueller 1996). The remnant masses can be de- WW15 (Figure 8). This indicates that our discussion above seemsto be quite valid, although we really need 1D results for WHW15to draw a more solid conclusion. mpacts of collective neutrino oscillations on supernova explosions 11
Table 3
2D simulationsModel Dimension R ν h ǫ ν x i L ν t s Explosion E ∞ diag M t = t s M ∞ [km] [MeV] [10 erg s − ] [ms] [10 erg] [ M ⊙ ] [ M ⊙ ]NH13A 2D — — — — Yes ∼ . Yes < . Yes < . Yes
Yes < . Yes termined by the mass ejection driven by the neutrinoheating, which range in 1.1-1.5 M ⊙ depending on the pro-genitors. For our 2D model of the 15 M ⊙ progenitor, wefound a set of the parameters that produces an explosionwith a canonical supernova energy close to 10 erg andat the same time leaves behind a remnant mass close to ∼ . M ⊙ . Our results suggest that collective neutrino os-cillations have the potential to solve the supernova prob-lem if they occurs. These effects should be explored ina more self-consistent manner in hydrodynamic simula-tions.Here it should be noted that the simulations in thispaper are only a very first step towards more realistic su-pernova modeling. For the neutrino transfer, we omittedthe cooling of heavy lepton neutrinos and the inelasticneutrino scattering by electrons. These omissions lead toan overestimation of the diagnostic energy and also theyshould relax the criteria for explosion. The ray-by-rayapproximation may lead to an overestimation of the di-rectional dependence of the neutrino anisotropies. A full-angle transport will give us a more correct answer (seeOtt et al. 2008; Brandt et al. 2011). Moreover, due tothe coordinate symmetry axis, the SASI develops prefer-entially along the axis; it could thus provide a more favor-able condition for the explosion. As several exploratorysimulations have been done recently (e.g., Iwakami et al.2008; Scheidegger et al. 2008; Nordhaus et al. 2010), 3Dsupernova models are indeed necessary also to pin downthe outcomes of the spectral swapping.Finally we briefly discuss whether the oscillation pa-rameters taken in this paper are really valid in viewsof recent work whose focus is on clarifying the still-veiled nature of collective neutrino oscillations. Follow-ing Duan et al. (2010), there are at least two conditionsfor the onset of collective neutrino oscillations in the caseof inverted neutrino mass hierarchy.The first criteria should be satisfied in the so-calledbipolar regime of the collective oscillation. In the regime,the neutrino number density should exceed the criticalvalue, n ¯ ν e , crit ≃ √ χ − ∆ m √ G F h ǫ ¯ ν e i≃ . × cm − (cid:18) . χ (cid:19) (cid:18)
15 MeV h ǫ ¯ ν e i (cid:19) , (12)where χ is the fractional excess of neutrinos over antineu-trinos, ∆ m is the characteristic mass-squared splitting(a typical value of ∼ . × − eV is employed here), χ Time after bounce [ms]
Figure 14.
Time evolution of χ = F ν e /F ¯ ν e −
1, where F ν i is thenumber flux of ν i . and G F is Fermi coupling constant. By using our simula-tion results, we can estimate χ which is often treated asa parameter (typically ∼ χ ≃ F ν e /F ¯ ν e − F ν x , where F ν i is the number flux of ν i . From Figure 14, it can be seenthat χ ∼ r ∼ n ¯ ν e = L ¯ ν e πr c h ǫ ¯ ν e i∼ . × cm − (cid:18) L ¯ ν e erg s − (cid:19) (cid:18)
100 km r (cid:19) × (cid:18)
15 MeV h ǫ ¯ ν e i (cid:19) , (13)therefore, the first condition is satisfied .The second criteria is related to the decoherence of col-lective oscillations by matter. In order to overwhelm thesuppression by the decohenrence, the following conditionshould be satisfied n ¯ ν e , crit ∼ n e , (14)where n e is the number density of electrons where the Even if χ is as small as χ ∼ .
01 due to the inclusion of ν x ,the criteria could be marginally satisfied. Y ¯ ν e , crit ∼ Y e . (15)In our 1D simulation, Y ¯ ν e ∼ (0 . − . × Y e for 100 km . r . r sh , where r sh is the shock radius . Since this con-dition is barely satisfied, the collective oscillations in re-ality could modify the spectrum to some extent betweenheavy-lepton neutrinos and electron/anti-electron neu-trinos, however the full swapping assumed in this studymay be exaggerated. Very recently , Chakraborty et al.(2011a,b) pointed out that the matter effect could fullysuppress the spectral swapping in the accretion phase us-ing 1D neutrino-radiation hydrodynamic simulation dataof Fischer et al. (2010). However, the current under-standing of the collective oscillation is not completedand calculations in this field employ several assump-tions (e.g., single angle approximation) (but see alsoDasgupta et al. 2011, for more recent work). To drawa robust conclusion, one needs a more detailed study including the collective neutrino flavor oscillation to thehydrodynamic simulations in a more self-consistent man-ner, which we are going to challenge as a sequel of thisstudy.We thank to K. Sumiyoshi, H. Suzuki, S. Yamada,T. Yoshida for stimulating discussions. Numerical com-putations were in part carried on XT4 at CfCA of theNational Astronomical Observatory of Japan. ML aresupported by the Swiss National Science Foundation un-der grant No. PP00P2-124879 and 200020-122287. Thisstudy was supported in part by the Japan Society forPromotion of Science (JSPS) Research Fellowships (YS),the Grants-in-Aid for the Scientific Research from theMinistry of Education, Science and Culture of Japan(Nos. 19104006, 19540309 and 20740150), and HPCIStrategic Program of Japanese MEXT. APPENDIX
CODE VALIDITY
Conservation of Energy
In this section, we demonstrate the conservation of physical quantities using the spherical collapse model (NH13).Figure 15 depicts the evolution of total binding energy by gravity (red line), total internal energy (green), total kineticenergy (blue), total trapped-neutrino energy (magenta), total energy leaked by neutrinos (cyan), and variation ofoverall energy (black dashed), respectively. The gravitational energy and total energy are negative and absolute valuesare shown. The gravitational energy and internal energy dominate (with different sign) and reach ∼ erg soonafter bounce. Despite such an enormous energy change, the total energy varies only within ∼ × erg so that theviolation of energy conservation remains < . ν e emission) that is consistent with the shock stall. We have monitored these values in a 2Dsimulation and obtained a similar level of energy conservation.
0 0.1 0.2 0.3 0.4 0.5 E n e r g i e s [ e r g ] Time after Bounse [s]Gravitational Binding EnergyInternal EnergyKinetic Energy Trapped-Neutrino EnergyLeaked Energy by NeutrinosVariation of Total Energy
Figure 15.
Time evolution of gravitational energy (red), internal energy (green line), kinetic energy (blue line), trapped-neutrino energy(magenta line), released energy by neutrinos (cyan line), and summation of these energies (black dashed line). These quantities aredetermined by integration with respect to volume included in our simulation except for released energy by neutrinos (magenta), which is R ( L ν e + L ¯ ν e ) dt . Since gravitational energy and total energy are negative, the absolute values are shown. The violation of total energy(dashed line) remains < × erg, which is ∼ .
03% of gravitational energy and internal energy after bounce ( ∼ erg). Outside the shock, Y ¯ ν e > Y e is achieved due to rapid densitydecrease. In fact they posted their papers on astro-ph after our submis- sion. mpacts of collective neutrino oscillations on supernova explosions 13
Comparison with AGILE
Here, we present the result of our numerical simulation in spherical symmetry and compare with the result of AGILE-IDSA code (Liebend¨orfer et al. 2009). AGILE (Adaptive Grid with Implicit Leap Extrapolation) is an implicit generalrelativistic hydrodynamics code that evolves the Einstein equations based on conservative finite differencing on anadaptive grid. We employ a one-dimensional version of our ZEUS-2D code that has been developed to performmultidimensional supernova simulations.We compare the evolution of a 13 M ⊙ star of Nomoto & Hashimoto (1988) in Newtonian gravity from precollapsemodel to 100 ms after bounce. We find good agreement between the results of the ZEUS-2D and AGILE during theearly postbounce phase when the neutrino burst is launched and the accretion shock expands to its maximum radius.The hydrodynamic quantities are shown in following figures. D e nd it y [ g c m - ] Mass Coordinate [M ⊙ ] ZEUS-2DAGILE 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 E n t r opy [ k B ] Mass Coordinate [M ⊙ ]ZEUS-2DAGILE -9e+09-8e+09-7e+09-6e+09-5e+09-4e+09-3e+09-2e+09-1e+09 0 1e+09 2e+09 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 V e l o c i o t y [ cc m s - ] Mass Coordinate [M ⊙ ] ZEUS-2DAGILE 0.25 0.3 0.35 0.4 0.45 0.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 Y e Mass Coordinate [M ⊙ ]ZEUS-2DAGILE Figure 16.
Density (left top), entropy (right top), velocity (left bottom), and electron fraction (right bottom) as a function of enclosedmass for the result of ZEUS-2D (red lines) and AGILE (green lines). The comparison is shown at the time just after the bounce. Adifference is seen in the entropy profile, which comes from the difference of shock capturing scheme. D e nd it y [ g c m - ] Mass Coordinate [M ⊙ ] ZEUS-2DAGILE 0 1 2 3 4 5 6 7 8 9 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 E n t r opy [ k B ] Mass Coordinate [M ⊙ ]ZEUS-2DAGILE -5e+09-4e+09-3e+09-2e+09-1e+09 0 1e+09 2e+09 3e+09 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 V e l o c i o t y [ cc m s - ] Mass Coordinate [M ⊙ ] ZEUS-2DAGILE 0.25 0.3 0.35 0.4 0.45 0.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 Y e Mass Coordinate [M ⊙ ]ZEUS-2DAGILE Figure 17.
Same as Fig. 16 but for the time at 1 ms after bounce. D e nd it y [ g c m - ] Mass Coordinate [M ⊙ ] ZEUS-2DAGILE 1 2 3 4 5 6 7 8 9 10 11 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 E n t r opy [ k B ] Mass Coordinate [M ⊙ ]ZEUS-2DAGILE -3e+09-2.5e+09-2e+09-1.5e+09-1e+09-5e+08 0 5e+08 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 V e l o c i o t y [ cc m s - ] Mass Coordinate [M ⊙ ]ZEUS-2DAGILE 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 Y e Mass Coordinate [M ⊙ ]ZEUS-2DAGILE Figure 18.
Same as Fig. 16 but for the time at 100 ms after bounce.REFERENCESBethe, H. A., & Wilson, J. R. 1985, ApJ, 295, 14Blondin, J. M., Mezzacappa, A., & DeMarino, C. 2003, ApJ, 584, 971Brandt, T. D., Burrows, A., Ott, C. D., & Livne, E. 2011, ApJ, 728, 8Bruenn, S. W., Mezzacappa, A., Hix, W. R., Blondin, J. M., Marronetti, P., Messer, O. E. B., Dirk, C. J., & Yoshida, S. 2009, inAmerican Institute of Physics Conference Series, Vol. 1111, American Institute of Physics Conference Series, ed. G. Giobbi,A. Tornambe, G. Raimondo, M. Limongi, L. A. Antonelli, N. Menci, & E. Brocato, 593–601Buras, R., Janka, H., Rampp, M., & Kifonidis, K. 2006, A&A, 457, 281Burrows, A., Dessart, L., Livne, E., Ott, C. D., & Murphy, J. 2007a, ApJ, 664, 416Burrows, A., Hayes, J., & Fryxell, B. A. 1995, ApJ, 450, 830Burrows, A., Livne, E., Dessart, L., Ott, C. D., & Murphy, J. 2007b, ApJ, 655, 416Chakraborty, S., Fischer, T., Mirizzi, A., Saviano, N., & Tomas, R. 2011a, arXiv:1104.4031Chakraborty, S., Fischer, T., Mirizzi, A., Saviano, N., & Tomas, R. 2011b, arXiv:1105.1130Colgate, S. A., & White, R. H. 1966, ApJ, 143, 626Dasgupta, B. 2010, arXiv:1005.2681Dasgupta, B., Dighe, A., & Mirizzi, A. 2008, Physical Review Letters, 101, 171801Dasgupta, B., O’Connor, E. P. & Ott, C. D. 2011, arXiv:1106.1167
Demorest, P. B., Pennucci, T., Ransom, S. M., Roberts, M. S. E., & Hessels, J. W. T. 2010, Nature, 467, 1081Duan, H., Fuller, G. M., & Carlson, J. 2008, Computational Science and Discovery, 1, 015007Duan, H., Fuller, G. M., & Qian, Y. 2010, Annual Review of Nuclear and Particle Science, 60, 569Esteban-Pretel, A., Pastor, S., Tom`as, R., Raffelt, G. G., & Sigl, G. 2007, Phys. Rev. D, 76, 125018Fern´andez, R. 2010, ApJ, 725, 1563Fischer, T., Whitehouse, S. C., Mezzacappa, A., Thielemann, F., & Liebend¨orfer, M. 2010, A&A, 517, A80+Fogli, G., Lisi, E., Marrone, A., & Mirizzi, A. 2007, Journal of Cosmology and Astro-Particle Physics, 12, 10Foglizzo, T., Galletti, P., Scheck, L., & Janka, H. 2007, ApJ, 654, 1006Fryer, C. L., & Warren, M. S. 2002, ApJ, 574, L65—. 2004, ApJ, 601, 391Guilet, J., Sato, J., & Foglizzo, T. 2010, ApJ, 713, 1350Herant, M., Benz, W., Hix, W. R., Fryer, C. L., & Colgate, S. A. 1994, ApJ, 435, 339Iwakami, W., Kotake, K., Ohnishi, N., Yamada, S., & Sawada, K. 2008, ApJ, 678, 1207Janka, H. 2001, A&A, 368, 527Janka, H., & Mueller, E. 1996, A&A, 306, 167Kitaura, F. S., Janka, H., & Hillebrandt, W. 2006, A&A, 450, 345Kotake, K., Sato, K., & Takahashi, K. 2006, Rep. Prog. Phys., 69, 971Kotake, K., Sawai, H., Yamada, S., & Sato, K. 2004, ApJ, 608, 391Langanke, K., et al. 2003, Physical Review Letters, 90, 241102Lattimer, J. M., & Prakash, M. 2007, Phys. Rep., 442, 109Lattimer, J. M., & Swesty, F. D. 1991, Nuclear Physics A, 535, 331Liebend¨orfer, M., Messer, O. E. B., Mezzacappa, A., Bruenn, S. W., Cardall, C. Y., & Thielemann, F. 2004, ApJS, 150, 263Liebend¨orfer, M., Mezzacappa, A., Thielemann, F.-K., Messer, O. E., Hix, W. R., & Bruenn, S. W. 2001, Phys. Rev. D, 63, 103004Liebend¨orfer, M., Whitehouse, S. C., & Fischer, T. 2009, ApJ, 698, 1174Marek, A., & Janka, H. 2009, ApJ, 694, 664Masada, Y., Takiwaki, T., & Kotake, K. 2011, submitted to ApJM¨uller, B., Janka, H., & Dimmelmeier, H. 2010, ApJS, 189, 104Murphy, J. W., & Burrows, A. 2008, ApJ, 688, 1159Nomoto, K., & Hashimoto, M. 1988, Phys. Rep., 163, 13Nordhaus, J., Burrows, A., Almgren, A., & Bell, J. 2010, ApJ, 720, 694Obergaulinger, M., Aloy, M. A., & M¨uller, E. 2006, A&A, 450, 1107Obergaulinger, M., & Janka, H.-T. 2011, arXiv:1101.1198O’Connor, E., & Ott, C. D. 2011, ApJ, 730, 70Ohnishi, N., Kotake, K., & Yamada, S. 2006, ApJ, 641, 1018Ott, C. D., Burrows, A., Dessart, L., & Livne, E. 2008, ApJ, 685, 1069Raffelt, G. G., & Smirnov, A. Y. 2007, Phys. Rev. D, 76, 081301Rampp, M., & Janka, H. 2002, A&A, 396, 361Sagert, I., Fischer, T., Hempel, M., Pagliara, G., Schaffner-Bielich, J., Mezzacappa, A., Thielemann, F., & Liebend¨orfer, M. 2009,Physical Review Letters, 102, 081101Scheck, L., Kifonidis, K., Janka, H., & M¨uller, E. 2006, A&A, 457, 963Scheidegger, S., Fischer, T., Whitehouse, S. C., & Liebend¨orfer, M. 2008, A&A, 490, 231Shlomo, S., Kolomietz, V. M., & Col`o, G. 2006, European Physical Journal A, 30, 23Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753Strack, P., & Burrows, A. 2005, Phys. Rev. D, 71, 093004Sumiyoshi, K., Yamada, S., Suzuki, H., Shen, H., Chiba, S., & Toki, H. 2005, ApJ, 629, 922Suwa, Y., Kotake, K., Takiwaki, T., Whitehouse, S. C., Liebend¨orfer, M., & Sato, K. 2010, PASJ, 62, L49Suwa, Y., Takiwaki, T., Kotake, K., & Sato, K. 2007a, ApJ, 665, L43—. 2007b, PASJ, 59, 771—. 2009, ApJ, 690, 913Suzuki, T. K., Sumiyoshi, K., & Yamada, S. 2008, ApJ, 678, 1200Takiwaki, T., Kotake, K., & Sato, K. 2009, ApJ, 691, 1360Thompson, T. A., Burrows, A., & Pinto, P. A. 2003, ApJ, 592, 434Thompson, T. A., Quataert, E., & Burrows, A. 2005, ApJ, 620, 861Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Reviews of Modern Physics, 74, 1015Woosley, S. E., & Weaver, T. A. 1995, ApJS, 101, 181Yakunin, K. N., et al. 2010, Classical and Quantum Gravity, 27, 194005Yamada, S. 2000, Phys. Rev. D, 62, 093026 < ε ν > L ν [ e r g s - M e V ] M [M ⊙⊙