Toward connecting core-collapse supernova theory with observations: I. Shock revival in a 15 Msun blue supergiant progenitor with SN 1987A energetics
DDraft version August 6, 2018
Preprint typeset using L A TEX style emulateapj v. 5/2/11
TOWARD CONNECTING CORE-COLLAPSE SUPERNOVA THEORY WITH OBSERVATIONS: I. SHOCKREVIVAL IN A 15 M (cid:12)
BLUE SUPERGIANT PROGENITOR WITH SN 1987A ENERGETICS
Timothy Handy , Tomasz Plewa , and Andrzej Odrzywo(cid:32)lek Draft version August 6, 2018
ABSTRACTWe study the evolution of the collapsing core of a 15 M (cid:12) blue supergiant supernova progenitor fromthe core bounce until 1.5 seconds later. We present a sample of hydrodynamic models parameterizedto match the explosion energetics of SN 1987A.We find the spatial model dimensionality to be an important contributing factor in the explosionprocess. Compared to two-dimensional simulations, our three-dimensional models require lower neu-trino luminosities to produce equally energetic explosions. We estimate that the convective engine inour models is 4% more efficient in three dimensions than in two dimensions. We propose that thegreater efficiency of the convective engine found in three-dimensional simulations might be due to thelarger surface-to-volume ratio of convective plumes, which aids in distributing energy deposited byneutrinos.We do not find evidence of the standing accretion shock instability nor turbulence being a keyfactor in powering the explosion in our models. Instead, the analysis of the energy transport in thepost-shock region reveals characteristics of penetrative convection. The explosion energy decreasesdramatically once the resolution is inadequate to capture the morphology of convection on large scales.This shows that the role of dimensionality is secondary to correctly accounting for the basic physicsof the explosion.We also analyze information provided by particle tracers embedded in the flow, and find thatthe unbound material has relatively long residency times in two-dimensional models, while in threedimensions a significant fraction of the explosion energy is carried by particles with relatively shortresidency times.
Subject headings: hydrodynamics — instabilities — shock waves — supernovae: general INTRODUCTION
There is substantial controversy over the importanceof physics processes involved in core-collapse supernova(ccSN) explosions (see, e.g., Janka et al. 2012; Burrows2013, and references therein). One problem is a multi-tude of relevant processes which include hydrodynamics,gravity, nuclear equation of state, and neutrino-matterinteractions. Despite much effort in the last four decades,no model accounting for these physics effects exists thatconsistently produces energetic ccSNe explosions for suit-able progenitor models. Early theories of ccSNe explo-sions such as neutrino-driven convection (see, e.g., Bur-rows & Goshy 1993; Janka & Mueller 1996; Mezzacappaet al. 1998) and the more recent standing accretion shockinstability (SASI) (Blondin et al. 2003) remain the cen-terpiece of computational core-collapse supernova stud-ies.The SASI mechanism has been discovered in the courseof numerical simulations in the work by Blondin et al.(2003). Subsequent SASI studies included theoreticalanalysis (Laming 2007; Foglizzo 2009) and numerical in-vestigations (Ohnishi et al. 2006; Blondin & Mezzacappa2007; Scheck et al. 2008). This shock instability is char-acterized by global, low order oscillations of the super-nova shock which generate substantial non-radial flow
Corresponding author: [email protected] Department of Scientific Computing, Florida State Univer-sity, Tallahassee, FL 32306, U.S.A. Marian Smoluchowski Institute of Physics, Jagiellonian Uni-versity, Reymonta 4, 30-059 Cracow, Poland in the post-shock region. These fluid motions result inincreased residence times of the accreted material andpossible increase the heating effect by neutrinos emittedby the nascent proto-neutron star.Since SASI is only one of the participating processes,it is often times difficult to unambiguously quantify itscontribution to the overall explosion process. Only re-cently has careful analysis been attempted to quantifythe role of SASI and compare its contribution to the othermajor participating process of neutrino-driven convec-tion. In particular, SASI may only be present in progen-itors with higher masses when neutrino-driven convectionmight be suppressed by increased accretion rates (M¨ulleret al. 2012, see also Bruenn et al. (2013)). Furthermore,the situation is additionally complicated by effects dueto model dimensionality observed in some simulations(Murphy & Burrows 2008; Nordhaus et al. 2010; Hankeet al. 2012; Dolence et al. 2013; Couch 2013; Takiwakiet al. 2013). Those findings are, however, not confirmedby all studies (Janka et al. 2012). (In passing, we wouldlike to note that is conceivable that the high frequencyof SASI observations in two-dimensional simulations canbe due to the assumed symmetry of such models. In ei-ther case, disentangling this kind of effects has not beenattempted in a systematic manner.) There is also littleknown about the difference in the efficiency of convec-tive heat transport between two and three dimensions,especially in the context of core-collapse supernovae. Al-though three-dimensional ccSNe explosion models havebeen available for some time (see, for example, Fryer & a r X i v : . [ a s t r o - ph . S R ] J a n Handy, Plewa, & Odrzywo(cid:32)lekWarren 2002), no systematic study aimed at comparingthe convection efficiency between two and three dimen-sions exists.In this work, we investigate the energetic explosions ofa 15 M (cid:12) post-collapse WPE15 model of Bruenn (1993) inone, two, and three spatial dimensions. This model hasbeen used extensively in the past by Kifonidis and his col-laborators (Kifonidis et al. 2003, 2006; Gawryszczak et al.2010, see also Wongwathanarat et al. (2013)), mainly inapplication to post-explosion mixing of nucleosyntheticproducts. Here our focus is solely on the explosion mech-anism. To this end, we tune the parameterized neutrinoluminosity such that the explosion energy near saturationmatches the observationally constrained energetics of SN1987A. In the process, we have constructed a databasecontaining more than 250 individual model realizations,differing in model dimensionality, parameterized neu-trino luminosity, and the small amplitude noise addedto the initial velocity. Obtaining such a large databaseof models was required given the extreme sensitivity ofmodels to perturbations, and one cannot draw conclu-sions about the nature of the explosion process exceptin a statistical manner. Furthermore, the observed dif-ferences between individual model realizations reflect thecombined model sensitivity and may provide informationabout the role of, and coupling between, participatingphysics processes.In Section 2 we describe the physical and numericalmodels adopted in our study. In Section 3 we presentoverall properties of the explosion models. In particu-lar, we focus on dimensionality effects, provide estimatesof neutron star kick velocities, and analyze the flow dy-namics in the context of neutrino-driven convection andSASI. Next, in Section 4, we examine the role of convec-tion in more detail. We analyze the energy flow in gainregion using energy flux decomposition, and introduceand apply new methods for characterizing the morphol-ogy of the post-shock flow. Section 5 discusses the po-tential role of turbulence in the explosion process, whilein Section 5.3 we use the Lagrangian view of the flowand track the energy histories of a large sample of fluidparcels in the gain region in order to characterize con-ditions leading up to shock revival. Finally, Section 6contains our discussion of results and conclusions of thiswork. METHODS
Hydrodynamics
We model the conservation of mass, momentum, en-ergy, and electron fraction as ∂ρ∂t + ∇ · ( ρ u ) = 0 , (1) ∂ρ u ∂t + ∇ · ( ρ u ⊗ u ) + ∇ P = − ρ ∇ Φ + Q M , (2) ∂ρE∂t + ∇· [ u ( ρE + P )] = − ρ u ·∇ Φ+ Q E + u · Q M , (3) ∂ρY e ∂t + ∇ · ( ρY e u ) = Q N . (4)where ρ is the mass density, u the fluid velocity, P thethermal pressure, Φ the gravitational potential, E = u / e the specific total energy (with e being the spe-cific internal energy), and Y e is the electron fraction. Thesource terms Q M , Q E , and Q N account for the effects ofneutrino-matter interaction (see, e.g., Janka & Mueller(1996)).Equations (1)–(4) are evolved using the neutrino-hydrodynamics code of (Janka & Mueller 1996; Kifoni-dis et al. 2003; Janka et al. 2003). The hydrodynamicssolver of this code is built upon the Piecewise ParabolicMethod (PPM) of Colella & Woodward (1984), and ap-plies Strang splitting (Strang 1968) to evolve multidi-mensional problems.The numerical fluxes inside grid-aligned shocks arecalculated using the HLLE Riemann solver of Einfeldt(1988) to avoid the problem of even-odd decoupling (Ki-fonidis et al. 2006); otherwise, the hydrodynamic fluxesare calculated using the Riemann solver of Colella & Glaz(1985) is used. Equation of State
We adopt the equation of state (EOS) outlined byJanka & Mueller (1996). This EOS consists of contri-butions from free nucleons, α -particles, and a heavy nu-cleus in nuclear statistical equilibrium. Nuclei and nucle-ons are treated as ideal, nonrelativistic Boltzmann gases,with electrons and positrons contributing as arbitrarilydegenerate and arbitrarily relativistic, ideal Fermi gases.Thermal effects between photons and massive particlesare taken into account. Nuclear statistical equilibriumwas assumed for T (cid:38) . ∼ . × K). Goodagreement has been reported between this EOS and theLattimer & Swesty (1991) EOS for ρ (cid:46) × g cm − (Janka & Mueller 1996). Gravity
The Newtonian gravitational potential field is treatedas a composite of a point mass field (for the excised PNS),and a distributed field (for the mass on the mesh). Thecomposite field is computed using a multipole expansionof the Poisson equation for gravity, ∇ Φ Newton = 4 πGρ, (5)where G is Newton’s gravitational constant. Sphericallysymmetric effects due to general relativity are accountedfor using the effective relativistic potential of Rampp &Janka (2002) such that,Φ final = Φ Newton − Φ DNewton + Φ DT OV , (6)where Φ DNewton and Φ DT OV are the spherically symmetricNewtonian and Tolman-Oppenheimer-Volkoff potentials,respectively. See Kifonidis et al. (2003); Marek et al.(2006); Murphy & Burrows (2008) for the details of thismethod. Thus, only in two-dimensions we account for de-viations from spherical symmetry in the Newtonian partof the potential. Otherwise, both the Newtonian andgeneral relativistic contributions are spherically symmet-ric.
Numerical grid
All simulations were performed in spherical geome-try. In the radial direction we utilize 450 logarithmicallyhock revival with SN 1987A energetics 3spaced zones. The inner boundary is located at a time-dependent location to mimic contraction of the collaps-ing proto-neutron star core (see Section 2.9 for details ofthis parameterization). The outer boundary is located at4 × cm. The choice of outer boundary radius ensuresthat the total mass on the mesh remains essentially con-stant throughout the simulation. In multidimensionalmodels we excise a 6 ◦ cone around the symmetry axisin order to minimize potential numerical artifacts asso-ciated with the symmetry axis (see Scheck (2007) forjustification of this approach).Multidimensional models have angular resolution of 2 ◦ and 3 ◦ in 2D and 3D, respectively. Additionally, a rep-resentative 3D model was selected to conduct a seriesof coarse resolution simulations at 6 ◦ , 12 ◦ , and 24 ◦ (cf.Section 6). Boundary conditions
We impose reflecting boundary conditions at the sym-metry axis. In three dimensions, we use periodic bound-aries in the azimuthal direction. The outer radial bound-ary is transmitting (zero gradient), while the innerboundary is an impenetrable wall with the radial pres-sure gradient matching the hydrostatic equilibrium con-dition.In order to account for continued contraction of theneutron star, we parameterize the inner boundary posi-tion, r ib , according to r ib ( t ) = r iib − exp ( − t/t ib )] (cid:16) r iib /r fib − (cid:17) , (7)as described in Janka & Mueller (1996). The contractionof the proto-neutron star is parameterized by a charac-teristic timescale, t ib . We make this parameterizationunder the assumption that even in fully consistent mod-els, uncertain physics such as the neutrino-matter inter-action rates and the (stiffness of) the nuclear equation ofstate will inevitably lead to uncertainties in the cool-ing of the neutron star, and therefore its contractiontimescale. Larger values of contraction timescale cor-respond to stiffer nuclear equations of state (slow con-traction), while smaller timescales lead to softer nuclearequations of state (fast contraction).In the present work, we consider families of both fastand slow contracting PNS core models, denoted as SCand FC models, respectively, to reflect the above un-certainty. The first family are the slowly contractingproto-neutron star models with r fib = 15 km and t ib = 1s. The second family are fast contracting proto-neutronstar models with r fib = 10 . t ib = 0 .
25 s. In bothcases we use r iib = 54 .
47 km.
Neutrinos
In our implementation of neutrino energy depositionwe use a light-bulb approximation (Janka & Mueller1996) with modifications introduced later by Scheck et al.(2006) (see below). Accordingly, we impose an isotropic,parameterized neutrino luminosity, L tot, ν . The initial to-tal neutrino luminosity is chosen such that it asymptoti-cally represents the gravitational binding energy releasedby the neutron star core (Janka & Mueller 1996; Scheck et al. 2006; Scheck 2007),∆ E ∞ core = (cid:90) ∞ L tot, ν h ( t ) dt = 3 L tot, ν t ib , (8)The neutrino luminosity varies in time as L ν x ( r ib , t ) = L tot, ν K ν x h ( t ) , (9)where ν x ≡ ν e , ν e , ν µ , ν µ , ν τ , ν τ . The K ν x terms rep-resent the constant, fractional contributions from eachneutrino type to the total neutrino luminosity. The h ( t )function describes the temporal evolution of the neutrinoluminosity and is given by Scheck et al. (2006) h ( t ) = (cid:40) t ≤ t ib ,( t ib /t ) / if t > t ib . (10) Tracer Particles
We utilize passive tracer particles to provide additionalanalysis of the flow structure. Particles are initially dis-tributed uniformly in Lagrangian mass coordinate. Eachparticle represents a fluid parcel with a mass of 1 × − M (cid:12) in two dimensions and 1 × − M (cid:12) in three dimen-sions (thus, the total number of particles is 7 . × and7 . × in two and three-dimensional models, respec-tively). We note that our choice of distributing particlesuniformly in Lagrangian mass results in large numbers ofparticles near the inner boundary due to high densitiesthere. Consequently, not all particles sample the mostdynamically active parts of the flow. Therefore, we limitour analysis of particle-based data to those inside thegain region. Initial conditions
For the initial conditions we use the post-collapsemodel WPE15 ls(180) of Bruenn (1993) based on the 15M (cid:12) blue supergiant progenitor model by Woosley (1988).We map the spherically symmetric collapse model to thegrid and in multidimensions we add small velocity per-turbations with relative amplitude 1 × − of the radialvelocity. Tuning explosion energies
In the current work we are interested in explosion mod-els with energetics (with energy per unit mass) matchingthat of SN 1987A. Estimates of the SN 1987A explosionenergetics are summarized in Table 1). Based onthis summary, we adopt
E/M = 0 . × erg M (cid:12)− ,which for the adopted progenitor model corresponds toan explosion energy of approximately 1 . × erg.We define the total model positive binding energy as E exp = (cid:88) ρ (cid:20) u + ε + Φ (cid:21) , (11)where the summation is taken over grid cells for whichthe binding energy is positive, 1 / u + ε +Φ >
0. Here u isthe velocity magnitude, ε is the internal energy, and Φ isthe gravitational potential. We define the explosion timeas the moment when the positive binding energy exceeds1 × erg (Janka & Mueller 1996). From that point intime on the total positive binding energy is identified asthe explosion energy. Handy, Plewa, & Odrzywo(cid:32)lek Table 1
Estimated explosion energetics of SN 1987A.E/M M ej E exp (10 erg/M (cid:12) ) (M (cid:12) ) (10 erg)Arnett (1987) · · · · · · · · · · · · · · · ± · · · ± ± · · · · · · Blinnikov (1999) 0.75 ± ± ± ± In order to match the model energetics with that ofSN 1987A, we computed a large number of trial explosionmodels with various neutrino luminosities recording theirexplosion energies at the final time, t = 1 . GENERAL PROPERTIES OF EXPLOSION MODELS
We obtained a large database of supernova explo-sion models in one-, two-, and three-dimensions. Thedatabase contained 26 spherically symmetric models, 1692D models, and 61 3D models. The individual modelrealizations differed in mesh resolution, parameterizedneutrino luminosity, and random perturbation pattern.Table 2 provides a summary of a subset of the databaseof explosion models that most closely match the en-ergetics of SN 1987A. The subset contains 10 models intwo-dimensions for both slow and fast contracting fami-lies, and 5 slow contracting models in three-dimensions.The model explosion times vary from about 200 msin the case of slow contracting models to slightly above100 ms in the case of fast contracting models. Everyfamily of models shows intrinsic dispersion in both theexplosion times and the explosion energies. For example,the explosion times for slow contracting models vary byabout 30% in 2D, and about 5% in 3D. The variationsin explosion times are comparatively modest in the caseof fast contracting models in 2D, with observed varia-tions not exceeding 10%. The corresponding variationsin the energetics are about 20% in 2D for slow and fastcontracting models. The observed variations of our 3Dmodels does not exceed 4%.We discuss the accretion rates and shock radii in Sec-tions 3.1 and 3.4, respectively.
Effects of neutrino luminosity on explosion energy
Figure 1 shows the explosion energy in models with
Figure 1.
Explosion energy for models tuned to the energeticsof SN 1987A (saturation at ≈ . × erg). Thick curves de-note the average over all models in the group; thin curves denotethe respective minimum and maximum values. Note that the fastcontracting proto-neutron star models begin unbinding material inroughly half the explosion time of the standard contracting models. both slow and fast contracting PNS cores. We show themean explosion energy and the explosion energy envelopefor 2D FC models (dash-dot), 2D SC models (dash), and3D SC models (solid). In all cases, the explosion energyrises rapidly at the onset of explosion and approximatelysaturates by the final time. In the case of SC models,the explosion energy threshold is reached between 140 msto 200 ms post-bounce depending on model realization.The two dimensional models tend to explode faster andreach their final energies later than the three dimensionalmodels. The FC models explode significantly earlier, atapproximately 110 ms on average. The short explosiontimes found in FC models, as compared to SC models,are not unexpected. This is because significantly greaterneutrino luminosities, approximately by a factor of two,are required to energize the explosion in the presence ofthe deeper gravitational potential well in the FC models.Regardless of the contraction times of the PNS core, bythe final time, t f = 1 . . × erg.It is instructive to discuss our models in terms of thecritical luminosity curve (Burrows & Goshy 1993). Weshall note that finding the curve was not the goal of ourhock revival with SN 1987A energetics 5 Table 2
Parameters and main properties of the explosion models.Model a L ν e b t exp c E exp d ˙ M exp e r exps (10 erg/s) (ms) (10 erg) (M (cid:12) /s) (km)Slow contracting models - 2DM194A 1.943 161 0.973 0.282 657M194B 1.943 184 1.038 0.265 765M194C 1.943 202 0.933 0.254 820M194D 1.943 186 1.025 0.260 847M194E 1.943 167 1.145 0.274 744M194F 1.943 200 0.925 0.258 755M194G 1.943 166 1.080 0.276 759M194H 1.943 207 0.937 0.248 867M194I 1.943 194 1.105 0.260 732M194J 1.943 148 1.054 0.288 721Slow contracting models - 3DM187A 1.871 189 1.066 0.259 814M187B 1.871 193 1.034 0.256 822M187C 1.871 184 1.054 0.265 793M187D 1.871 187 1.036 0.260 818M187E 1.871 182 1.058 0.264 805Fast contracting models - 2DM352A 3.528 109 1.074 0.315 697M352B 3.528 114 1.037 0.311 711M352C 3.528 115 1.089 0.310 724M352D 3.528 112 0.995 0.312 707M352E 3.528 110 0.990 0.311 702M352F 3.528 114 1.139 0.313 707M352G 3.528 110 0.968 0.316 667M352H 3.528 107 1.221 0.317 666M352I 3.528 117 0.962 0.305 772M352J 3.528 113 1.019 0.309 729 a Parameterized electron neutrino luminosity. The model anti-electron neutrino luminosity is 7.5%larger. b Explosion time determined when the total positive binding energy is greater than 1 × erg. c Total positive binding energy (i.e. explosion energy) evaluated at t = 1 . d Mass flow rate evaluated immediately upstream from the shock at t = t exp . e Average shock radius evaluated at t = t exp . study, and therefore, we are not in a position to providethe exact relation. Our results only provide an upperlimit on the critical neutrino luminosity for the (relativelysmall) range of accretion rates present in our models (seeTable 2).Figure 2 provides a comparison of our explosion mod-els with the results of critical neutrino luminosity stud-ies conducted by other groups. In the figure, the re-sults of Nordhaus et al. (2010), Hanke et al. (2012), andCouch (2013) are shown with open circles, triangles, andsquares, respectively, while our data are shown with solidcircles; the model dimensionality is indicated with linestyle with solid, dashed, and dash-dotted lines connect-ing models obtained in 1D, 2D, and 3D, respectively.In the left panel, the neutrino luminosity is shown as afunction of the accretion rate, which is measured by in-tegrating the mass flux through the surface located justupstream of the shock at the onset of explosion.This figure indicates that our choice to constrain thefinal explosion energy by observations required a rela-tively narrow range of neutrino luminosities, comparedto other studies. We also find that the required neutrinoluminosity systematically decreases as the dimensional-ity of the model increases. Neutrino luminosities foundin our models only provide an upper limit for the criti- cal luminosity for the mass accretion rates present in oursimulations. Furthermore, explosion times found in ourmodels are shorter than those of other groups for similarneutrino luminosities. The simplest explanation for thisfinding is that the neutrino heating is more efficient in ourmodels. However, we cannot exclude the possibility thatthis observed, systematic difference is due to differencesin the progenitors. Detailed investigation of possible dif-ferences between our findings and those of other groupsis beyond the scope of the current work. Neutron star recoil
Since we excise the neutron star from the grid, a spe-cial method has to be used to estimate its recoil velocity.To this end, we use the approach of Scheck et al. (2006),and exploit momentum conservation. In this method,the momentum of the neutron star balances the momen-tum of the surrounding envelope (total momentum onthe mesh in our simulations). An additional correctioncan be applied to the neutron star recoil velocities byconsidering that the neutron star will continue accretingmaterial beyond the final time in our simulations. Weestimate the corrected velocity as v ∞ H = v . H + v − . H , (12) Handy, Plewa, & Odrzywo(cid:32)lek Figure 2.
Comparison of our explosion models with the results of critical neutrino luminosity studies conducted by other groups. Theresults of individual studies are shown with different symbols, while model dimensionality is represented by line style (1D: solid; 2D:dashed; 3D: dash-dotted). (left panel) The dependence of the critical neutrino luminosity on the accretion rate at the time of explosion.We observe a systematic decrease in the required neutrino luminosity as the model dimensionality increases. (right panel) Parameterizedelectron neutrino luminosity as a function of explosion time. For a given dimensionality, our models tend to explode sooner and requirelower neutrino luminosities. where v . H is the recoil velocity obtained in our simula-tions at t f = 1 . v − . H is the approximate cor-rection due to late time accretion onto the neutron star.To calculate the velocity correction, we assumed the rateof the recoil velocity change obtained by Scheck et al.(2006) in his B18 model series between the times t = 1 .
5s and t = 3 . (cid:28) v − . S v . − S (cid:29) ≈ . . (13)Then the recoil velocity correction has been computed as v − . H = (cid:28) v − . S v . − S (cid:29) v . − H , (14)In the above formulae superscripts denote post-bounceexplosion times, and the subscripts S and H denote thevalues from Scheck et al. (2006) and the current work,respectively.Figure 3 shows the evolution of recoil velocities in ourmodels until the final time (left panel), and the estimatedsaturation velocities based on Equation 12 (right panel).Our results indicate relatively modest recoil velocitiescompared to Scheck et al. (2006). The approximate finalrecoil velocities are in the range between 100 km/s and300 km/s, with the lowest and highest recoil velocitiesof 30 km/s and 475 km/s. These results are in qualita-tive agreement with the results reported by Scheck et al.(2006, Figure 20) and Wongwathanarat et al. (2013, B-series models in Table 2). Gain region characteristics
In this section we report the results following the anal-ysis of a quasi-steady state period in the pre-explosionepoch. If no such evolutionary stage is reached then theentire pre-explosion phase could be considered a tran-sient phenomenon with little importance for the ener-getic core-collapse supernova explosions considered here.If the opposite is true however, we can adopt methods
Figure 3.
Neutron star kick velocity for our SC models, shownwith dashed and solid lines, respectively. The outset shows the esti-mated saturation velocities based on estimates utilizing the resultsof Scheck et al. (2006). appropriate for analysis of quasi-steady state flows andapply them to analyze the dynamics of the gain region.We identify the quasi-steady state period in the post-shock flow evolution by calculating the mass containedin the gain region, M gain . Its time-dependence is shownin the left panel of Figure 4 for our set of multidimen-sional models. As one can see, the mass in the gainregion for the case of SC models (shown with solid lines)changes only slightly between t = 100 ms and t = 150ms. A similar period of modest increase in the mass ofthe gain region can be found in the case of FC models(shown with dashed lines) for times between t = 30 msand t = 80 ms. We note that in both cases, the momentwhen the mass in the gain region increases precedes theexplosion times, which is consistent with a scenario inwhich the shock revival process occurs over an extendedperiod of time in a quasi-steady state fashion. Also notehock revival with SN 1987A energetics 7 Figure 4. (left panel) Mass in the gain region. The evolution of mass in the gain region is shown with solid lines for SC models (thin solidlines in 2D and thick solid lines in 3D) and dashed lines for FC models. Note that the mass in the gain region is, on average, greater inFC models than in SC models. Also, after the initial transient oscillations, the mass in the gain region stabilizes, indicating that evolutionin the gain region has reached a quasi-steady state period. (right panel) Shock accretion rate. The line types are associated with modelfamily and model dimensionality as in the left panel.Note that after the initial period of fast accretion the accretion rates progressivelydecrease in both families of models. The accretion rates are, on average, greater in SC models than in FC models after the initial transient.See Section 3.3 for details. the presence of strong oscillations in the gain region massat early times in both families of models. Since, as wediscuss below (Section 3.4), the shock radius is increasingsteadily at these early times, the observed oscillations inthe gain region mass are either due to density changes ofthe material entering the gain region or the changes inthe position of the gain radius during these times. Thoseoscillations, however, have a transient character and donot play any role in the subsequent evolution of the sys-tem, and cease as soon as the flow in the gain regionbecomes multidimensional.The mass accretion rate (shown in the right panel ofFigure 4) as measured at the position of the shock frontshows a rapid decline during the first 50 ms of the sim-ulation time after which it steadily decreases at a pro-gressively slower rate. The evolution of accretion ratedoes not differ significantly between SC and FC models.This is expected, as we consider only a single progeni-tor model. At later times when the quasi-steady state isreached, the accretion rates are between 0 . − .
35 M (cid:12) s − for SC models and between 0 . − .
67 M (cid:12) s − forFC models. The FC models are characterized by consis-tently lower accretion rates than the SC models for timeslater than 50 ms.To better understand the dynamics in the gain regionin the context of neutrino heating, we consider the ad-vection timescale, τ adv , τ adv = M gain ˙ M , (15)where M gain is the mass in the gain region, and ˙ M isthe accretion rate. The advection time scale is the char-acteristic time a fluid parcel spends in the gain region.This quantity neglects multidimensional effects, whichare known to be important in the process of shock revival.In the same context, we consider the heating efficiency, η , η = τ adv (cid:82) ρ ˙ QdV , (16)
Figure 5.
Advective timescale for the gain region. The oscil-latory behavior of the advective timescale at early times for SCmodels likely has the same origin as the transient oscillations ob-served in the gain region mass. Note that the advective timesare similar between model realizations in the same family. Thedifferences in advective times for different model realizations areparticularly large for 2D SC models that develop a few large-scalenonuniformities in the gain region. On average, the advective timesare longer for FC models than for SC models. See Section 3.3 fordetails. where ˙ Q is the net neutrino energy deposition rate, andthe integral in the above equation is taken over the grav-itationally bound ( u + ε + Φ <
0) material inside thegain region. The quantity in the denominator of theabove equation is the heating timescale, i.e. the time ittakes to heat a gravitationally bound fluid parcel so thatit becomes unbound. The heating efficiency is a measureof the competition between the advection and heatingprocesses. In particular, η > ≈
100 km) than in our models ( ≈ t = 20 ms and t = 50 ms in FC models, and t = 50ms and t = 90 ms in SC models. This increase in theheating efficiency is associated with the onset of convec-tion prior to the quasi-steady state phase (see Section 4).At later times, after the quasi-steady state is established,the heating efficiency continues to steadily increase in theSC models while the increase in the heating efficiency issignificantly greater in the FC models. Nevertheless, theheating efficiency continues to evolve in a qualitativelysimilar manner in both families of models. As we notedbefore, the similarities between SC and FC models canbe easily explained by the fact that we consider only asingle progenitor model in this work.On the quantitative level, the heating efficiency in theFC models is greater than one from the beginning of theevolution. On the other hand, the heating becomes effi-cient in the SC models only after the quasi-steady state isestablished. We conclude that a strong, immediate neu-trino heating is required in order to produce energeticexplosions. Shock evolution
The shock revival process can also be studied by an-alyzing the time evolution of the shock radius. In par-
Figure 6.
Heating efficiency for the gain region. The evolutionof heating efficiency is shown for SC and FC models with solid anddashed lines, respectively. The results from the 2D SC models areshown with thin solid lines. Note that the heating efficiency inthe FC models is always greater than one. This indicates that astrong, immediate neutrino heating is required in order to produceenergetic explosions in this case. On the other hand, the heat-ing becomes efficient in the SC models only once the quasi-steadystate is established. Again, as in the case of the mass in the gainregion and advective times, the evolution of heating efficiency isqualitatively similar between SC and FC families. ticular, prolonged periods of shock stagnation may indi-cate the onset of the standing accretion shock instability(SASI). Conversely, a steady increase in the shock radiusmay indicate that other processes (e.g. neutrino-drivenconvection) operate efficiently. Additionally, the unbind-ing of shocked matter is directly dependent on radius(through the gravitational potential), and large shockradii may aid in reaching the explosion threshold. Fi-nally, hydrodynamic perturbations present at early timescan be imprinted on the shock and affect the morphol-ogy of the supernova ejecta. To this end, in the followingdiscussion we use the shock aspect ratio to quantify theexplosion asymmetry.The evolution of the average shock radius, R s , in ourmodels is shown in Figure 7 with solid and dashed linesfor the SC and FC models, respectively. The shock radiusincreases steadily in both families of models, although theshock expansion is about 25% faster in the FC models.At the time when the early convection sets in, the shockexpands to about 270 km in the FC models and about200 km in the SC models. By the time the shock isrevived, its average radius for different model realizationsvaries between 645 km and 855 km in the 2D SC models,and between 790 km and 810 km in the 3D SC models.We believe the greater variation in the average shockradius in the 2D models compared to 3D may hint atdifferences in the dynamics of the gain region between2D and 3D. Moreover, 2D FC models display a similaramount of variation in the average shock position as 2DSC models. Also, in that case, the average shock positionat the time of explosion varies between 640 km and 730km. Overall, the average shock radius at the time ofexplosion is similar in SC and FC families. One shouldkeep in mind, however, that the FC models explode, onaverage, twice faster than the SC models (c.f. Table 2).The shock aspect ratio, r maxs /r mins , found in our mod-hock revival with SN 1987A energetics 9 Figure 7.
Evolution of the average shock radius. The averageshock radius is shown with solid and dashed lines for SC and FCmodels, respectively, with 2D SC models shown with thin solidlines. Shaded areas denote regions where ( t exp , r s ) pairs for SCand FC model families are located. Note that the FC and SCmodels explode at similar shock radii but at significantly differenttimes. The shock expansion rate for FC models is, on average,twice that of the SC models. Figure 8.
Shock aspect ratio for our models. Thick lines denotethe average over all models in the group. Thin lines denote theminimum and maximum values for the particular families. Ourresults show relatively mild asymmetries ( r maxs /r mins ≈ . els (see Figure 8) appears relatively modest compared tothe shock asymmetry reported in other studies (Hankeet al. 2012; M¨uller et al. 2012). Although an extremeshock aspect ratio of 2 was found in 2D SC simulations,the average aspect ratios are only 1.25 in 2D and 1.4 in3D. This should be compared to the shock aspect ratiosbetween 3 and 4 we estimated based on the data pre-sented by M¨uller et al. (2012).The most striking feature of the evolution of the shockaspect ratios in our models are strong asymmetries ob-served in the 2D SC models at early times. These asym-metries seem to rapidly decrease for times >
400 ms. Incontrast, the shock aspect ratios in the 2D FC modelsappear to evolve more smoothly and the most deformedFC models seem to retain their asymmetric character astime progresses. Qualitatively, however, the degree to which the shock is deformed in 2D models is modest andon average is ≈ .
3. In 3D models, the average degreeof shock asymmetry at late times is somewhat greater( ≈ . l = 1 mode developsin those models (cf. Figure 10(a)) and results in a strongasymmetry. In other 2D SC models, the higher ordermodes make the shock more spherical (cf. Figure 10(b)).The situation is similar in the case of 3D SC models,in which the post-shock region shows richer convectivestructure and evolves on a similar timescale as in 2D.The evolution of shock radius also provides informa-tion that is helpful in the context of the standing accre-tion shock instability (SASI). Specifically, SASI mani-fests itself as low-order ( l = 1 ,
2) oscillations of the shockradius (Foglizzo et al. 2006; Laming 2007; Scheck et al.2008; Foglizzo 2009). Therefore, the first step in the SASIanalysis of the shock revival is to decompose the shockradius in terms of spherical harmonics. The expansioncoefficients are given by a lm = (cid:90) Ω Y ∗ lm ( θ, φ ) f ( θ, φ ) d Ω , (17)where the spherical harmonic Y ∗ lm is normalized by a fac-tor (cid:113) l +14 π ( l − m )!( l + m )! , and f ( θ, φ ) = r s ( θ, φ ).In order to enable direct comparison between the 2D( m = 0) and 3D ( m = − l . . . l ) cases, we consider a suit-ably normalized contribution of the m -mode coefficients(Ott et al. 2013), α l = 1 a (cid:118)(cid:117)(cid:117)(cid:116) l (cid:88) m = − l a lm . (18)The evolution of α and α coefficients is shown in Figure9. The larger values of the coefficients indicate that the2D SC models, on average, exhibit stronger variabilitythan their three-dimensional counterparts. However, wealso observe a certain dichotomy among 2D SC modelswith one subgroup showing variability distinctly largerthan the remaining set of models. 2D SC models alsoshow relatively weak l = 1 mode contributions. pertur-bations, while others rise to a moderate fraction of theaverage shock position. The 2D FC models (data notshown in Figure 9) show a similar degree of variabilityas the 2D SC models at all times.Prior to explosion, the 2D models exhibit at most oneor two weak oscillations. This is a quantitatively differentbehavior than found in the simulations of M¨uller et al.(2012), who reported strong, multiple cycles leading upto shock revival. In addition, the shock perturbationsin our models are weak compared to those consideredas evidence for SASI. Furthermore, we do not find qual-itative differences between the 2D and 3D models, aswe mentioned above. Note that most evidence for SASIpresented in the literature is essentially restricted to 2D,non-exploding models.0 Handy, Plewa, & Odrzywo(cid:32)lek Figure 9.
Evolution of the leading coefficients in the spherical harmonic decomposition of the shock radius. Only the results of thedecomposition for SC models are shown. The runs of the α (left panel) and α (right panel) coefficients are presented with solid anddashed lines for 3D and 2D models, respectively. Prior to the onset of convection ( t ≈
50 ms), the shock front remains essentially sphericallysymmetric. The emergence of low order perturbations around t ≈
80 ms is due to buoyant convective plumes deforming select sections ofthe shock front rather than SASI. Note that the degree of shock perturbation is relatively modest compared to some recent results (see,e.g., M¨uller et al. 2012).
In our SC models, the low order modes emerge at about t = 75 ms, shortly before the quasi-steady state is estab-lished in the gain region. Visual inspection of the flowmorphology during the quasi-steady state provides noevidence for large-scale “sloshing” motions in the gainregion, considered a defining signature of SASI. We con-clude it is unlikely the SASI plays any important role inthe evolution of our models. Morphology of the gain region
As soon as the neutrino-driven convection sets in andthe related perturbations reach the shock, one faces a dif-ficult problem of disentangling various physics processesparticipating in the shock revival, including fluid insta-bilities and neutrino-matter interactions. Analysis of theoverall flow morphology is the first step in the processof understanding the explosion mechanism. It providesinitial evidence for the physics of fluid flow participatingin re-energizing the stalled shock (Herant et al. 1994). Italso may provide evidence for the possible role of modelparameters, such as assumed symmetries and discretiza-tion errors, e.g. near the symmetry axis (Scheck et al.2006; Gawryszczak et al. 2010), on the flow dynamics.In this section we present the morphological evolution ofour SC models in two and three dimensions. In partic-ular, we are interested in identifying when the fluid flowinstabilities imprint their structure on the just-formedinner core of the supernova ejecta.Entropy pseudocolor maps of the gain region for twotwo-dimensional and two three-dimensional SC modelsat their respective explosion times are shown in Figure10. Prior to reaching the explosion threshold, the 2DSC models evolve to have either dipolar ( l = 1, Figure10(a)) or quadrupolar ( l = 2, Figure 10(b)) m = 0 ejectamorphology. These flow features emerge when strongdownflows in the gain region are formed shortly after theonset of convection. Furthermore, our models indicatethat these flow structures do not evolve once they set inand continue to persist after the explosion is launched,evolving in a self-similar fashion at later times (see be-low). Additionally, these downflows carry accreted ma- terial deep into the gain region close to the gain radius.This behavior does not seem to be as extreme in 3D SCmodels, because a comparable amount of accreted ma-terial is transported through many more downflows (seeFigure 10(c) and (d)). The presence of many downflowsin three-dimensions reflects the fact that the flow struc-ture is inherently different in 2D and 3D. We quantifythis behavior later in Section 4.The entropy distribution in the 3D SC model M187Ais shown in Figure 11 for select times. At early times,the flow evolves away from initially spherically symmet-ric accretion, with the first signs of convective instabilityemerging around t = 50 ms (Figure 11(a)). At this time,the shock is not affected by convection. At intermediatetimes ( t = 100 ms, Figure 11(b)), the initial convectiveplumes gradually merge into larger structures and be-gin deforming large segments of the shock front. As timeprogresses, the process of bubble merging continues (Fig-ure 11(c)) and the shock eventually is launched around t = 189 ms (Figure 11(d)).The large-scale morphology of the post-shock regionappears essentially frozen after the explosion commences.For example, three ejecta plumes seen in the southernhemisphere (Figure 11(d)), can still be easily identifiedmore than 300 ms later (Figure 11(e)). However, dur-ing the same period the convective bubbles in the north-ern hemisphere show significant evolution and appear tomerge into a single structure. This morphology persistsuntil the end of the simulation and the central region be-comes filled with the neutrino-driven wind (Figure 11(f)). CONVECTION
Our analysis of the explosion process presented aboverevealed that the explosion is preceded by a relativelyshort (compared to other studies) quasi-steady state in-side the gain region (Section 3.3) during which the shockonly slowly expands and SASI does not seem to a ma-jor role (Section 3.4). We also found that the heatingbecomes efficient shortly before the quasi-steady state isestablished (Section 3.3).There exists a large body of evidence from numericalhock revival with SN 1987A energetics 11 (a)
M194Jt exp = 148 ms r [ k m ] z [ km ] (b) M194Ct exp = 202 ms r [ k m ] z [ km ] (c) M187At exp = 189 ms y [ k m ] x [ km ] (d) M187Bt exp = 193 ms y [ k m ] x [ km ] Figure 10.
Entropy distribution in the post-shock region at the explosion time for select two-dimensional and three-dimensional models.Entropy distribution is shown with pseudocolor maps for models M194J (panel a), M194C (panel b), M187A (panel c), and M187B (panel d).The entropy distribution maps for two-dimensional models (panels (a), (b)) are reflected across the symmetry axis. For three-dimensionalmodels, the entropy is shown for a slice through the computational domain at the equatorial plane. In two dimensions models tend to haveeither dipolar (panel a) or quadrupolar (panel b) ejecta structures which result in an oblate shape of the shock. In three dimensions theejecta appear to have more regular structure, and the shock is less deformed. See Section 3.5 for discussion. simulations of core-collapse supernova explosions indicat-ing that multidimensionality, and specifically the pres-ence of convection, decreases the requirements for theneutrino luminosity produced by the contracting core ofthe proto-neutron star (Janka & Mueller 1996; Nordhauset al. 2010; Hanke et al. 2012; Couch 2013). In thissection we focus on characterization and quantificationof the effects of neutrino-driven convection during thequasi-steady state of the gain region. We introduce andapply novel methods dedicated to investigating this as-pect of the supernova explosion process.
Minimum resolution requirements
We perform a series of three-dimensional simulationsin order to assess the dependence of model physics op-erating inside the gain region on numerical resolution.In particular, we are interested in finding conditions re-quired to suppress neutrino-driven convection on largescales. If indeed convection is a critical component driv-ing the explosion, one could expect that by suppressingconvection (by whatever means) the explosion would notoccur, even in models with neutrino luminosities thatin other situations are sufficiently large enough to re-vive the shock. This expectation is strongly supportedby numerous calculations which show that, for example,much higher luminosities are required to produce an ex-plosion in one-dimension than in multidimensions (see,2 Handy, Plewa, & Odrzywo(cid:32)lek (a)t = 50 ms y (cid:2) c m (cid:3) x (cid:2) cm (cid:3) (b)t = 100 ms y (cid:2) c m (cid:3) x (cid:2) cm (cid:3) (c)t = 150 ms y (cid:2) c m (cid:3) x (cid:2) cm (cid:3) (d)t = 189 ms y (cid:2) c m (cid:3) x (cid:2) cm (cid:3) (e)t = 500 ms y (cid:2) c m (cid:3) x (cid:2) cm (cid:3) (f)t = 1500 ms y (cid:2) c m (cid:3) x (cid:2) cm (cid:3) Figure 11.
Distribution of entropy in the 3D SC model M187A. The entropy is shown with pseudocolor maps for a slice through thecomputational domain at the equatorial plane. The entropy distribution is shown at t = 50 , , , , , and 1500 ms in panels (a)-(f).The development of convection can be seen as early as t = 50 ms (panel a). The flow structure becomes progressively more complex atintermediate times (panels (b) and (c)) and the explosion is launched at t = 189 ms (panel (d)). Relatively little change in the structureof the gain region takes place at later times (panels (e) and (f)). Table 3
Characteristics of the M187A modelseries with varying angular resolution.Angular t exp E exp ( t = 250 ms)Resolution (ms) (10 erg)2 ◦
183 0.1143 ◦
183 0.1076 ◦
175 0.10612 ◦
198 0.07724 ◦
216 0.028 e.g., Janka & Mueller (1996), Section 3.3, and Figure 2for discussion and compilation of relevant recent results).In our series of 3D simulations, we keep the radial res-olution fixed while the angular resolution is graduallydecreasing by a factor of 2 between the models, from 3 ◦ (our base resolution) down to 24 ◦ . In addition, we com-pute a 2 ◦ resolution model. For each model, we recordedthe explosion times and explosion energies at t = 250ms. The basic characteristics of models obtained in thisseries is given in Table 3. The inspection of datashown in the table reveals that the explosion times andthe explosion energies change little as long as the angu-lar resolution is not worse than 6 ◦ . For coarser angularmeshes, we observe a significant increase in the explo-sion times and decrease in the explosion energies. Noexplosions were found at resolutions coarser than 24 ◦ .It is conceivable that the observed changes in globalmodel characteristics should be correlated with changesin the flow structure inside the gain region. Figure 12 shows the entropy maps in a subset of the models usedin this resolution study. One can easily note profoundstructural changes in the morphology of the gain regionas the angular model resolution decreases. Large- andsmall-scale convective bubbles can clearly be seen in the3 ◦ resolution model (panel (a)). There are significantlyfewer small, convective bubbles present in the 6 ◦ model(panel (b)). In those two models the overall shock radiusappears comparable. In contrast, no well-defined, small-scale bubbles can be identified in the 12 ◦ or 24 ◦ models(shown in panels (c) and (d), respectively). In addition,the shock radius is visibly smaller in the 24 ◦ model.We interpret the observed dependence of energeticsand explosion timing on angular resolution describedabove as follows. Consider that the hydrodynamic solverused in our simulations is the PPM scheme (Colella &Woodward 1984). PPM is nominally third order (andpractically second order) accurate in space, and usespiece-wise parabolic interpolation to describe profiles ofthe hydrodynamic state. Therefore, the minimal numberof mesh resolution elements required to resolve a convec-tive bubble is about 3 mesh cells. In our 12 ◦ and 24 ◦ resolution models there are only 13 and 7 cells in angle,respectively. This implies our hydrodynamic solver canonly represent between 2 to 4 large-scale structures inthose models. The number of those structures does notrepresent the actual number of bubbles, as some meshcells are also required to describe the flow structure be-tween the bubbles. This explains the lack of bubbles onsmaller scales and the overall degraded appearance of theneutrino-driven convection in these cases.hock revival with SN 1987A energetics 13 (a)3 ◦ y [ k m ] x [ km ] (b)6 ◦ y [ k m ] x [ km ] (c)12 ◦ y [ k m ] x [ km ] (d)24 ◦ y [ k m ] x [ km ] Figure 12.
The distribution of entropy in the M187A series of models with varying angular resolution. The entropy is shown for aslice taken at the equatorial plane at t = 100 ms in models with angular resolutions of (a) 3 ◦ , (b) 6 ◦ , (c) 12 ◦ , and (d) 24 ◦ . Note theoverall appearance of convection is degraded as the angular resolution decreases, indicating a drop in the efficiency of the neutrino-drivenconvective engine. The shock radius is significantly smaller in the coarsest model. See Section 4.1 for discussion. We demonstrated through this resolution study thatthe neutrino-driven convection is a critical componentof the explosion mechanism in our models. Specifically,we have shown that one can turn an energetic, multi-dimensional model into a failed multidimensional modelsimply by degrading its angular resolution. Therefore,to correctly capture the efficiency of the convective en-gine one needs to resolve its basic components. This pic-ture is also consistent with the analysis of Herant et al.(1994), who argued that, independent of dimensionality,large scales will play the dominant role in convection.Furthermore, it is important to resolve not only struc-tures on the largest scales, but also the structures on ∼ ◦ are important. We conclude with the somewhatobvious statement that capturing the relevant physics ofthe problem requires adequate numerical resolution. Itis conceivable that the next generation of ccSNe models with much higher resolution will begin uncovering newphysics effects that cannot be observed in the currentgeneration of models due to their insufficient quality. Dynamics
One of the quantities of interest in the context of dy-namics of convection is the amount of mass containedinside rising, convective bubbles. To characterize thatquantity, we introduce the upflow mass fraction,ˆ m up ( t ) = (cid:82)(cid:82)(cid:82) gain ρ | u r > dVM gain , (19)where the integration is performed over the gain region,and takes into account only the fluid elements (gridzones) with positive radial velocity. Thus, the upflowmass fraction characterizes the mass inside the gain that4 Handy, Plewa, & Odrzywo(cid:32)lek Figure 13.
The temporal evolution of the upflowing mass in thegain region. The fraction of the total mass inside the gain regionwith ( u r >
0) is shown with dash-dotted, dashed, and solid linesfor 2D FC, 2D SC, and 3D SC models, respectively. See Section4.2 for discussion. is moved away from the proto-neutron star toward theshock. The evolution of the upflow mass fraction in ourmodels is shown in Figure 13. As one can see, the ini-tial behavior of the upflowing mass is characterized bya transient period until t ≈
50 ms. Soon after that, theamount of upflowing mass begins to steadily increase. Itlevels off during the quasi-steady state phase lasting ap-proximately between 85–140 ms for SC models (shownwith solid lines in Figure 13) and 65–105 ms for FCmodels (shown with dashed lines in Figure 13), and risesagain after the shock is launched. Interestingly, the netmass flux inside the gain region is close to zero in theSC models (the horizontal dashed line at ˆ m up = 0 . Figure 14.
The evolution of relative mass fluctuations insidethe gain region shortly before and during the quasi-steady statephase in the 3D SC explosion models. Fluctuations are calculatedrelative to a linear fit of the upflowing mass in the gain regionfrom the beginning of the quasi-steady state phase until the explo-sion time, where individual curves terminate. (thin solid) M157A;(medium solid) M157B; (thick solid) M157C; (dashed) M157D;(dotted) M157E. See Section 4.2 for discussion. outflowing mass in these models is non-monotonic. Onecan identify several short (∆ t ≈
10 ms) episodes duringwhich the amount of mass in upflows fluctuates. Thesefluctuations have an amplitude of several percent. Giventhat the typical number of convective bubbles, and there-fore also the number of downflows, in those models dur-ing the quasi-steady state phase is about 10 to 25 (seeSection 4.3), it is conceivable that perhaps individualdownflows are responsible for the observed fluctuations.
Characterization of convective structures
Having in mind the continuing debate regarding therole of dimensionality in numerical models of ccSNe ex-plosions, we seek information about the dependence oflarge-scale structure of convection inside the gain regionon model dimension. As the first step in our analysis, wewish to measure the net radial momentum over a fractionof the gain region, I mom ( θ, φ ) = r mid ( θ,φ )+ αw (cid:82) r mid ( θ,φ ) − αw ρu r dr >
00, otherwise (20)where r mid = ( r s + r g ) / r s and r g , re-spectively), w = r s − r g is the angle-dependent gain re-gion width, and α is a fractional offset between 0 and 0 . I mom = 1) ordownflows ( I mom = 0). Recall that our models are essen-tially free of possible SASI contributions, and thereforethis method should be applicable to any convective flows.Apart from the situation when the shock and gain ra-dius are constant, the radial extent of the integrationbounds in the above equation will vary with angle. In thisway, we avoid possible ambiguities that may arise nearthose two regions, and make the results of this methodhock revival with SN 1987A energetics 15 Figure 15.
Estimate of the total solid angle spanned by buoyant,rising material. The onset of convection at t ≈
50 ms results ina rapid increase in the spanned solid angle as the first convectiveplumes begin to form. During the quasi-steady convective phasethe plumes continue to grow. By the time the explosion is launchedthey cover roughly 2/3 of the total solid angle. See Section 4.3 fordiscussion. insensitive to flow asymmetries naturally developing es-pecially along the shock front. In the above equation weheuristically set α = 0 . t <
50 ms), the whole gain region oscillates radially. Asconvective structures begin to emerge, the solid angle oc-cupied by bubbles starts to increase. This process startsto slow down at the beginning of the quasi-steady statephase ( t ≈
80 ms). However, and unlike in the case ofother diagnostics we have discussed above, the evolutionof the solid angle occupied by bubbles does not provideany clear signature of the explosion time. Rather, thegrowth continues at a similar pace until about 150 msafter the explosion commences. It starts to level off soonafter that time, and by 400 ms after bounce the wholegain region is filled-in with upflowing material.Apart from large-amplitude oscillations, the evolutionof the solid angle occupied by bubbles in 2D SC models isqualitatively similar to that seen in 3D, with the averagevalue similar between the models. However, in a few 2Dcases persistent downflows develop (see also Section 3.5).For these models the solid angle occupied by bubblesdoes not reach 1 by the end of the period we analyzedhere. Note that our diagnostics can be use to identifysituations when strong accretion onto the proto-neutronstar continues at late times.Figure 16 shows the time evolution of the number ofbubbles inside the gain region averaged separately for 2Dand 3D SC models (shown with dashed and solid lines,
Figure 16.
Evolution of the number of rising bubbles inside thegain region. The number of bubbles is averaged separately for 2Dand 3D SC model realizations (shown with dashed and solid lines,respectively). For the clarity of presentation, the vertical scale islimited to (cid:104) N (cid:105) = 25. The maximum average number of bubblesfound in 3D is ≈ respectively). Initially, our method identifies only a sin-gle bubble or finds no bubble during the initial transient( t <
50 ms), in agreement with the radial oscillationsof the gain region during that period as we discussedearlier. Soon after convection sets in, the number of dis-tinct upflows sharply increases, reaching ≈
12 in 2D and ≈
200 in 3D. (Note, for the purpose of presentation wehave limited the scale in Figure 16 to 25.) Given thatat those early times the solid angle occupied by bubblesis small, we conclude that the bubbles are initially smalland grow in angular extent over time. By the time of theexplosion, the bubbles are roughly 4 times smaller in 3Dthan in 2D.Recall that in 2D the bubbles are truly tori rather thanquasi-spherical plumes. We recognize this may have cer-tain implications for the bubble dynamics (see, e.g., Ki-fonidis et al. 2003; Hammer et al. 2010; Couch 2013).Furthermore, the observed evolution of bubble sizes fromsmall to large (see also Herant et al. 1994) offers an in-teresting parallel between the evolution of the structureof neutrino-driven convection (and perhaps convection ingeneral), and that of bubble merger observed in a multi-mode Rayleigh-Taylor instability (see, e.g., Alon et al.1994; Miles 2004). The study of the process of mergingconvective bubbles is beyond the scope of the presentwork.The growth of convective bubbles after t = 50 ms co-incides with time at which the heating efficiency startsto increase (cf. Figure 6) and accumulation of mass inthe gain region (cf. Figure 4). We believe this indicatesthat early, small-scale convection begins to interact withthe incoming accretion flow, resulting in rather rapid in-crease in the gain region mass. Subsequently, the advec-tion time scale also increases. We believe this providesevidence for the direct dependence of the heating effi-ciency on the intensity of convection. Convective energy transport
In our discussion of convective energy transport, weadopt the approach originally introduced for the anal-ysis of stellar convection by Hurlburt et al. (1986). In6 Handy, Plewa, & Odrzywo(cid:32)lekthis approach, the individual components of the energytransport equation, Equation 3, are averaged over lat-eral directions. Then, one computes radial distributionsof deviations away from the lateral averages, f (cid:48) ≡ f − ¯ f .In the discussion here, we use the notation of Moc´aket al. (2009). We define the convective flux, F C , kineticflux, F K , buoyant work, P A , and expansion work, P P , asfollows: F C = (cid:73) u r ρ (cid:18) ε + Pρ (cid:19) (cid:48) r d Ω , (21) F K = (cid:73) u r ρ (cid:18) u · u (cid:19) (cid:48) r d Ω , (22) P A = − (cid:73) u r ρ (cid:48) ∂ Φ ∂r r d Ω , (23) P P = (cid:73) ( ∇ · u ) P (cid:48) r d Ω , (24)where the integrals are taken over the suitable solid angle(in our models with the excised 12 ◦ cone around the sym-metry axis this amounts to ≈ .
5% of the full solid an-gle). These terms describe the total energy transportedper unit time through a surface of a sphere of radius r .Figure 17 shows the time averages of convective andkinetic energy fluxes (left panel) and the buoyant and ex-pansion work terms (right panel) for the SC models. Themain characteristics of the energy transport follow fromthe direction of the energy transport and the relative con-tributions of individual terms. We begin our analysis atthe gain radius located at ≈
115 km in both 2D and 3Dmodels (marked with the solid, vertical line). There, theconvective flux is initially negative but rapidly increasesand reaches the maximum at ≈
15% of the radial ex-tent of the gain region. Further out, the convective fluxgradually decreases and ultimately ceases well upstreamof the shock. On the other hand, the kinetic energy fluxis negative below the gain radius and remains negativethrough most of the gain region. The observed run ofthe convective flux provides clear evidence for convectivetransport in which material heated near the gain radiusbuoyantly rises toward the shock.The relation between buoyancy and convection alsomanifests itself through the buoyant work, P A (shownwith gray lines in the right panel of Figure 17). Thework done by buoyancy becomes positive starting at thegain radius, indicative of buoyancy driving the convec-tion upwards. It first rapidly increases through the bot-tom layers of the gain region, and reaches its maximumat ≈
10% of the gain region extent. Further out, thework done by buoyancy gradually decreases and eventu-ally vanishes as expected at the shock, where the flowbecomes, on average, spherically symmetric.As we noted, the buoyant work is positive throughoutthe gain region. This can be caused by either overdensematerial sinking or underdense material rising. Sincethe work due to buoyancy is negative below the gainregion, therefore the opposite scenario with either over-dense, cold material is rising or underdense, hot materialis sinking. Since we do not observe sinking of the under-dense, hot material in that region, we conclude that thenegative work due to buoyancy indicates the presence ofrising overdense and cold matter. In our models, this rising, overdense material originates from the downflowsthat are turned around by the pressure gradient. Thisprocess is known as buoyant braking (Brummell et al.2002).The most conspicuous difference between 2D and 3Dlies in the distribution and value of the expansion worksource term, P P (shown with solid lines in the right panelof Figure 17). When integrated over time and space, thisterm represents the P dV work of fluid elements. Belowthe gain radius this term is negative, indicative of thework done in the process of compressing the fluid. Itchanges sign across the gain radius, and remains positivethrough the lower one-third of the gain region. This in-dicates that the fluid is expanding in that region. Wequantify the amount of work due to expansion by inte-grating Equation 24 over the portion of domain between r ≈
60 km and r ≈
250 km, where the expansion worksource term is significant and differs most between 2Dand 3D. Assuming a characteristic timescale of 100 ms,we estimate the (kinetic) energy change due to
P dV workto be ≈ × erg and ≈ − × erg in 2D and 3D,respectively, or about a few percent of the explosion en-ergy.We believe the difference in the P dV work, as we de-scribed above, reflects the fact that the structural in-tegrity of upflows and downflows in 2D is much greaterthan in 3D. The structural integrity of flow features de-pends on the surface-to-volume ratio. For a given vol-ume, 2D structures have less surface area than their 3Dcounterparts. Thus, flow structures in 2D are less sus-ceptible to (surface) perturbations. We alluded to thisproperty in our discussion of morphological differencesbetween 2D and 3D simulations in Section 3.5. Conse-quently, in 3D, the downflows do not penetrate as deep,which can be clearly seen in the left panel of Figure17, which shows that both the kinetic fluxes are nega-tive and of comparable magnitude in both 2D and 3Daround r ≈
250 km. However, as one moves into deeperlayers the kinetic flux much more rapidly decreases in2D than in 3D, and vanishes in the layers closer to theproto-neutron star surface. The convective fluxes in 2Dand 3D behave in a qualitatively similar way, althoughthey are positive and equal closer to the gain radius, at r ≈
150 km. From there inward, the convective fluxdecreases at a significantly greater rate across the gainradius and ceases closer to the proto-neutron star surfacein 2D than in 3D. Furthermore, the greater structural in-tegrity of upflows and their coherent nature in 2D makesthe expansion process dramatically more efficient thanin 3D. This is evidenced by much greater work done byexpansion above the gain radius, as we discussed earlier.As one can see in the right panel of Figure 17, thebuoyant and expansion terms are both negative belowthe gain region. This is because below the gain regionthe system is no longer energy conserving due to intenseneutrino heating. This situation is different, however,across the shock front. There, both terms are of differ-ent sign, which is consistent with the fact that the neu-trino heating is weak and the energy is approximatelyconserved in that region.Finally, we would like to note an interesting possibilitythat the difference in efficiency of the convective enginebetween 2D and 3D might be another consequence of thehock revival with SN 1987A energetics 17
Figure 17.
Decomposition of the total energy flux and the distribution of the related work source terms in the SC models. Quantitiesare averaged over time over the period between t = 110 ms and t = 120 ms. The individual curves represent averages taken over all modelrealizations for a particular model dimension (shown with dashed and solid lines for 2D and 3D models, respectively). The vertical linesdenote the approximate locations of the gain radius and the minimum shock radii for the two families. (left panel) Convective flux, F C ,and kinetic flux, F K . (right panel) Buoyant ( P A ) and expansion ( P P ) work. See Section 4.4 for discussion. difference in surface-to-volume ratio between those mod-els. This is because in 3D the surface area of the interfacebetween hot and cold material where the mixing, andthus also heat exchange, takes place is relatively larger.Quantifying this possibly important effect is, however,beyond the scope of this paper. TURBULENCE
Motivated by the recent discussion of the possible roleof turbulence in the process of shock revival (Murphy& Meakin 2011; Hanke et al. 2012), we present spectraof the kinetic energy in a lateral direction and turbulentReynolds stresses characteristic of the gain region for ourSC models. In particular, we will determine if the kineticenergy spectra are consistent with model dimensionality,and discuss the qualitative differences in the Reynoldsstresses found between 2D and 3D.
Spectra
We investigate the spectral properties of turbulenceby considering the spherical harmonic decomposition ofEquation 17 with f ( θ, φ ) = √ ρu θ , where u θ is the po-lar component of the velocity vector. Our choice of thisparticular velocity component hinges on the assumptionthat the flow is isotropic in the lateral directions, whichis supported by our analysis of the Reynolds stresses pre-sented below in Section 5.2.To ensure that the energy contained in the sphericalharmonic modes adequately represents the lateral kineticenergy, we include the factor √ ρ (Endeve et al. 2012).Following Hanke et al. (2012), we define the lateral ki-netic energy spectrum, E ( l ), as E ( l ) = l (cid:88) m = − l (cid:20)(cid:90) Ω Y ∗ lm ( θ, φ ) √ ρu θ d Ω (cid:21) , (25)where the integral is taken over the solid angle included inour simulations. In order to account for spatial and tem-poral variations of the lateral kinetic energy spectrum,we average the spectrum over a shell of thickness 25 kmcentered at the midpoint of the gain region ( r mid ≈ Figure 18.
Lateral kinetic energy spectra for SC models. Thespectra are shown with dashed and solid lines for 2D and 3D mod-els, respectively. The data are averaged spatially over a sphericalshell 25 km thick centered at 290 km, and over the time periodbetween t = 120 ms and t = 130 ms. Piecewise powerlaw functionsare fit separately to spectra for 2D and 3D model families, andare shown with dashed and solid lines for 2D and 3D, respectively.The fit process allows for optimal locations of the break points inthe piecewise powerlaw functions. See Section 5.1 for discussion. km). In addition, in our analysis the spectra are averagedover snapshots from the time interval between t = 120 msand t = 130 ms (soon after the quasi-steady state phase).In passing we note that there exists other choices of phys-ical quantities for the analysis of turbulence, for exam-ple density fluctuations (Borriello et al. 2013), which aremore suited for the analysis of neutrino-matter interac-tions rather than fluid flow dynamics.Figure 18 shows the lateral kinetic energy spectra forour SC models with 2D and 3D models shown withdashed and solid lines, respectively, averaged over timeand space as described above. We find that one canidentify three distinct spectral regions in the lateral ki-netic energy spectrum. The first region occupies largescales ( l (cid:46) a few), the intermediate region extends up8 Handy, Plewa, & Odrzywo(cid:32)lekto l (cid:46) several 10s, and the third region extends towardstill smaller scales ( l (cid:38) E ( l ) ∝ l α . In particular, it is expected that inthe intermediate regime, in which the energy is trans-ported from large to small scale (and possibly also in theopposite direction in 2D situations), a cascade developswith the power law exponent, α , chiefly dependent onthe particular characteristics of the system.To quantify the shape of the power spectra found inour simulations, we fit power law functions inside eachregion. We found that the power law exponent that bestcharacterizes the shape of the spectrum at large scalesis α ≈ − . α ≈ − . l ≈
35 with α ≈ −
4. At still smaller scales ( l (cid:38) (cid:46) l (cid:46) α = − . α = − . α = − /
3; Monin & I’Aglom (1971)).On a final note, the Reynolds number estimated for ourmodels is quite small ( Re (cid:39) . . . ≈ Turbulent Reynolds stresses
Figure 19 shows the time- and model realization-averaged turbulent Reynolds stress components for theSC models. The data shown in the figure were generatedas follows. First, the tensor components were computedfor individual model realizations using R ij = (cid:104) ρu i u (cid:48) j (cid:105) ρ , (26)where (cid:104)·(cid:105) denotes the averaging operation over solidangle, u (cid:48) j = u j − (cid:104) u j (cid:105) is the perturbation away fromthe background of the j -th component of velocity, and ρ = (cid:104) ρ (cid:105) is the background density. After averaging overthe solid angle, the stress components were then averagedover the period between t = 110 ms and t = 130 ms. Inthe final step, we used these averages to compute the Figure 19.
Distributions of Reynolds stresses in the gain regionfor the SC models. The data shown in the figure is obtained byfirst averaging stress distributions for individual model realizationsover the period between t = 110 ms and t = 130 ms, and then bycalculating mean values based on the time-averages separately for2D and 3D. The results for 2D and 3D families are shown with thinand thick lines, respectively. The dashed vertical lines at r ≈ model-averaged values of turbulent Reynolds stresses.In 3D, the radial (shown with the thick solid line in Fig-ure 19) and lateral stress components (shown with thickdash-dotted and dotted lines; their sum is shown with thethick dashed line), increase from the inner boundary andacross the gain radius (marked with a vertical, dashedline at r ≈
115 km). The radial component reaches themaximum around r ≈
300 km, roughly in the middle ofthe gain region, while the overall contribution of lateralstresses peaks slightly deeper at r ≈
255 km. The radialstress is about a factor of 2 greater than the sum of lat-eral stresses. Farther out, the stresses decrease and reachapproximately similar magnitude just below the shock.In the figure, the three dashed, vertical lines located be-tween 400 km and 600 km mark the innermost, average,and outermost shock radii. The stresses show a mild in-crease across the region occupied by the shock and even-tually vanish upstream of the shock inside the accretionflow. The sum of lateral stresses exceeds the radial stressby ≈
30% in a relatively narrow region behind the shock(at r ≈
430 km). Before we begin discussing the qual-itative and quantitative between the average turbulentReynolds stresses in 2D and 3D, we note that the radialdependence of stresses in 2D (shown with thin solid anddash-dotted lines for the radial and lateral stresses, re-spectively, in Figure 19) is quite similar to that in 3D.Moreover, flow appears roughly isotropic below the gainradius, with radial and lateral stresses in equipartition.We note that the stresses in 2D are much larger inthe region near the gain radius than in 3D. This canbe understood by recalling that the convection is mostvigorous in that region in our two-dimensional models(see Section 4.4). Farther out, the dominant, buoyantly-generated radial stresses undergo redistribution in lateraldirections (see, e.g., Murphy & Meakin (2011)). Thishock revival with SN 1987A energetics 19process of redistribution of the radial stress is expectedto be isotropic in the lateral directions, as observed inour 3D simulations (thick dash-dotted and dotted lines inFigure 19). By volume-integrating the radial averages ofthe stress components inside the gain region, we find thatin both 2D and 3D the integrated radial stress is largerby about 60% than the lateral stresses. Furthermore, theintegrated turbulent stresses are greater by about 30%in 3D than in 2D. The latter result may appear at firstsurprising, judging by the run of turbulent stresses shownin Figure 19. However it is important to recognize thatthe curves shown in the figure represent angular averagesrather than volume-integrated contributions. Therefore,even though the stresses are greater in 2D than in 3Dnear the gain radius, the larger stresses in 3D in the upperlayers of the gain region result in their overall greatervolume-integrated value.The fact that the volume-integrated turbulent stressesare greater in 3D than in 2D is consistent with thesurface-to-volume difference between 2D and 3D simu-lations (the argument we made in Section 4.4). In thisscenario, 3D structures are more susceptible to pertur-bations and the flow becomes disorganized over largerregions (in our case, a greater fraction of the radial ex-tent of the gain region). Specifically, full-width at half-maximum measure of the radial stress distribution is 25%greater in 3D than in 2D. (The corresponding volume fac-tor is still larger due to the 3D distribution being centeredat a higher radius.) Finally, we note that the process ofredistribution of the radial stress is appears more grad-ual in 3D than in 2D (compare the two curves connectedwith the double arrow in Figure 19). In fact, this mayonly be a misimpression. Note that the contribution oflateral stresses in 2D differs qualitatively from that in3D only in the lower half of the gain region. We findpreliminary evidence that this may be due to the differ-ence in the work done by expansion in that region, wherethe expansion work source term is large in 2D, while itis nearly zero in 3D (see Figure 17). We speculate thatthe expansion increases the velocity fluctuations and thuscontributes to the increase of the turbulent stresses.
Lagrangian analysis of the explosion mechanism
We believe that potentially important insights intothe core-collapse supernova explosion mechanism can begained through analyzing the history of fluid parcels asthey enter the gain region and participate in the revivalof the shock. Our motivation partially follows from therealization that the amount of evidence provided by theEulerian analysis appears insufficient to disentangle andclearly identify the role various physics plays in the ex-plosion process. To this end, we seeded our simulationswith a large number of tracer particles from the begin-ning of the simulations. Recall from Section 2 than inour models, each tracer particle represents a fluid ele-ment with a mass of 1 × − M (cid:12) in two-dimensions and1 × − M (cid:12) in three-dimensions. Residence times in the gain region
In our analysis we assume a particle resides inside thegain region if its position is between the shock radiusand the gain radius (both of which are angle- and time-dependent). We define the residency time, t res , as the Figure 20.
Distributions of tracer particles in the gain region asa function of particle residency time for SC models. This distribu-tions are shown with dashed and solid lines for 2D and 3D families,respectively, at t = 50 ms, 100 ms, and 200 ms. The data are aver-aged over all model realizations for a particular model dimension.For legibility, we smooth the curves with a boxcar method using atime window of 5 ms. Note that the scale differs for 3D models.See Section 5.3.1 for discussion. total time a tracer particle spends in the gain region.Thus, the maximum residency time of any particle neverexceeds the model time elapsed since the beginning ofthe simulation. We construct a distribution of particleresidency times by binning the residency times with aresolution of 1 ms.Figure 20 shows the model realization-averaged distri-butions of particle numbers for our SC models at t = 50ms, 100 ms, and 200 ms, as a function of particle resi-dency time. In the figure, the data for 2D and 3D modelsare shown with dashed and solid lines, respectively. Thetotal number of particles in the gain region are approx-imately 5 . × , 6 . × , and 7 . × in 2D and5 . × , 5 . × , and 7 . × in 3D, for progres-sively later evolutionary times. We note that particlesoccupying bins near t res = 0 should be interpreted asmaterial that has been just accreted (entered the gainregion), while the farther to the right the particle bin isthe earlier it entered the gain region.There is little difference in distributions of particlenumbers between the 2D and 3D models (shown withthick dashed and solid lines, respectively, in Figure 20)for particles that enter the gain region prior to the onsetof convection. This is expected, given that at those earlytimes the shock is nearly spherically symmetric and thepost-shock flow is essentially radial. These initial dis-tributions are expected to evolve due to both accretionof new material (particles passing through the shock),and loss of particles already residing in the gain regionthrough the gain radius and their settling onto the proto-neutron star. Therefore, one may expect certain differ-ences in the distributions of particle numbers at latertimes. First, notice that the number of particles withmaximum residency times is always smaller in 3D thanin 2D, except at the early time ( t = 50 ms). This meansthat a portion of particles with the longest residencytimes that are accreted through the gain radius is largerin 3D than in 2D. Second, this trend reverts at t ≈ t = 200 ms there are more particles with shortresidency times present in 3D than in 2D. One possibleexplanation for this behavior might be a systematic dif-ferences in accretion rates and shock radii between 2Dand 3D, but we found no evidence to support this possi-bility (cf. Section 3.3).Additional information about the origin and evolu-tion of particle number distributions can be obtained byexamining the spatial dependence of particle residencytimes. These distributions reveal that, in 3D, there ex-ists a population of particles with low residency timeslocated near the gain radius. The systematic differencesbetween particle number distributions in 2D and 3D cantherefore be explained by particles leaving and reenteringthe gain region across the gain radius in 3D. This leads tothe apparent deficit of particles with long residency timesand causes excess of particles with low residency times inthis case. We quantified the magnitude of this effect byexamining the contours of cumulative particle numbersat a given residency time. We found that between 10%and 20% of particles located near the gain radius partici-pate in this process. This estimate is consistent with thedifference between the particle number distributions seenin Figure 20. The above effect signifies a difference in theflow dynamics near the gain radius between 2D and 3D.In particular, we believe this behavior is due to a largeramount of mass (by 10–20%) involved in convection atthe bottom of the gain region, where the convection isdriven, and may be the reason for the greater efficiencyof convection in 3D (see Section 5.3.3 for discussion). Thermodynamic evolution of the shocked material
One quantity of particular interest in the context ofthe ccSNe explosion mechanism is the boost to the in-ternal energy of the material residing in the gain regionprovided by neutrino heating. This process can be ana-lyzed by considering the increase in the internal energy ofshocked fluid parcels (represented in our study by tracerparticles) as a function of the particle residency time.Figure 21 shows the model realization-averaged inter-nal energy gain per unit mass, (cid:104) ∆ e (cid:105) , as a function of par-ticle residency time for the 2D and 3D SC models (shownwith dashed and solid lines, respectively). We define theinternal energy gain of a tracer particle as the differencebetween the particle internal energy at the evolutionarytime shown and its internal energy when it enters thegain region. We find that at early times ( t (cid:46)
50 ms)there is essentially no difference in the way particles gainenergy in 2D and 3D (thick dashed and solid lines). Atlater times, we begin to observe modest at first ( t = 100ms) changes with the newly accreted material progres-sively gaining energy faster in 3D than in 2D (the solid,medium thick line representing 3D data lies always to theleft of the dashed, medium thick line representing the 2Ddata).In addition to the apparent asymmetry in the 3D dis-tribution at late times, we also find qualitatively newunique features, in particular at t = 200 ms. First, weidentify a region with almost constant increase in theenergy gain for particles that entered the gain region be-tween t = 50 ms and t = 100 ms. During this time, andas we discussed earlier in Section 5.3.1, the gain radiusshrinks faster in 3D than in 2D. This allows the particlesresiding close to the gain radius to gain comparatively Figure 21.
Distribution of the specific internal energy change fortracer particles in the gain region as a function of particle resi-dency time for SC models. The distributions are shown at t = 50,100, and 200 ms with thick, medium, and thin lines, respectively.The data are averaged over all model realizations for a particulardimension, and are shown separately with dashed and solid linesfor 2D and 3D. For legibility, we smoothed the curves with a box-car method using a time window of 5 ms. See Section 5.3.2 fordiscussion. more energy in 3D. This process is responsible for theobserved shape of the internal energy gain distributionfor residency times between t res ≈
100 ms and t res ≈ t = 200 ms(thin solid line) at t res ≈
100 ms. We believe this signi-fies the increase in the thermodynamic efficiency of theconvective engine in 3D. Third, there is a smaller humpfor a population of particles that entered the gain regionshortly before ( t res <
50 ms) the shock was launched.Finally, we would like to point out that the similarityin the shapes of the 3D distributions at late times isquite striking. The arrows in the figure connect specificfeatures in the energy gain distributions that are wellpreserved in 3D between t = 100 ms and t = 200 ms.This indicates that the particles which entered the post-shock region during the first 100 ms of the simulationgained internal energy at a similar rate during the processof shock revival. Unbinding of the shocked material
The explosion process necessarily requires gravitation-ally unbinding a substantial fraction of the collapsingcore. Alas, very little information about this process canbe found in the literature, perhaps with the exceptionof discussions of global quantities such as the explosionenergy. However, armed with the tracer particles we arenow in a position to investigate this process in more de-tail. In particular, the results discussed in the previoussection indicate that profound differences in the heatingefficiency exist between 2D and 3D models. Althoughwe are able to identify populations of particles respon-sible for those differences, our focus was more on theinterplay between the advection and neutrino heating ofhock revival with SN 1987A energetics 21
Figure 22.
Distribution of the positive specific gravitationalbinding energy for tracer particles in the gain region as a functionof particle residency time for SC models. The data are taken at t = 200 ms, and are averaged over all model realizations for a par-ticular model dimension. The results are shown with dashed andsolid lines for 2D and 3D, respectively. For legibility, we smooththe curves using a boxcar method with a time window of 5 ms.Note that in 2D, unbound material is composed primarily of mat-ter which entered the gain region prior to the quasi-steady phase.However, in 3D, there are substantial contributions from materialwhich entered the gain region during the period of vigorous con-vection. See Section 5.3.3 for discussion. the shocked material. In the following discussion of theprocess of gravitational unbinding of the shocked, stellarcore we also take into account the flow dynamics insidethe gain region.Figure 22 shows the model realization-average of thepositive part of the specific gravitational binding energy, E b , as a function of residency time for 2D and 3D fam-ilies of models at t = 200 ms (marked by dashed andsolid lines, respectively). The data indicates dramaticdifferences in the energetics of fluid elements in the gainregion between 2D and 3D. In 2D, the unbound portionof the gain region is almost exclusively made of materialthat entered the gain region during the first 100 ms of theevolution. However, in 3D one can identify two differentpopulations of particles. The first 3D population consistsof particles that entered the gain during the first 100 ms.Although this is the same material that is responsible forthe explosion in 2D, in 3D this material is substantiallyless energetic. A careful examination of the data shownin Figure 22 reveals that, indeed, around the explosiontime the explosion energies in 2D are greater than in 3D.This is consistent with the 2D models exploding, on av-erage, faster than 3D models. Also, as discussed in theprevious section, the particles that enter the gain regionearly on gain energy more efficiently in 2D than in 3D.It is only at later times when the particles in 3D show adramatic increase in their internal energies.In the scenario presented above, the explosion islaunched in 2D early on by a relatively stronger cen-tral neutrino source. In 3D, however, the more efficientconvective engine is able to launch the shock at neutrinoluminosities lower by ≈ SUMMARY & FUTURE WORK
We have presented an analysis of two- and three-dimensional models of early core-collapse supernova de-velopment in the collapsing core of a 15 M (cid:12) progenitor.In our study, the simulations start at shock bounce, con-tinue through the onset of neutrino-driven convection,quasi-steady state of the gain region, and early post-revival shock expansion until the energy of the explosionhad approximately saturated at 1.5 seconds. Our modelswere tuned to match the estimated explosion energeticsof SN 1987A by careful selection of the parameterizedneutrino luminosity. We considered the cases of slowand fast contraction rates of the proto-neutron star toreflect on uncertainty in the nuclear equation of state.We introduced and applied new diagnostic methods forthe explosion process, morphology and structure of con-vection, and energetics of material inside the gain region.Our conclusions are formulated based on a largedatabase of explosion models. This computationally ex-pensive approach is necessary since the core-collapse su-pernova problem involves highly nonlinear, strongly cou-pled physics which results in the extreme sensitivity ofrealistic computational models to perturbations. There-fore, conclusions of core-collapse supernova studies canonly be understood in a statistical sense. Moreover, suchensemble-based conclusions are potentially affected byresolution effects, especially in three-dimensions, due tothe high computational cost of individual model realiza-tions.The main findings of our work can be summarized asfollows: • We found that the same explosion energy can beobtained in three-dimensions with less energy thanin two-dimensions. This trend also holds betweentwo-dimensions and one-dimension. This result isin agreement with some previous studies. • The stochastic character of the explosion processresults in non-uniqueness of the relation betweenneutrino luminosity and explosion energy. In par-ticular, in the process of tuning the neutrino lumi-nosities, we were able to obtain equally energetic3D models with the neutrino luminosities differingby as much as 3%. • Dimensionality also contributes to the non-2 Handy, Plewa, & Odrzywo(cid:32)lekuniqueness of the above relation. In particular,equally energetic explosions can be obtained in 3Dusing neutrino luminosities smaller by 4% than in2D. This implies the efficiency of the convective en-gine is greater by that amount in 3D. • We proposed that the observed dependence of theefficiency of the convective engine on the model di-mensionality might be a consequence of the differ-ence in the surface-to-volume ratio of flow struc-tures between 2D and 3D. In this scenario, thegreater efficiency of the convective engine in 3Dis simply due to a larger surface-to-volume ratio inthat case, which helps the process of heat exchangebetween hot and cold flow structures. Quantifyingthis possibly important effect is, however, beyondthe scope of this paper. • We note that the relations between neutrino lumi-nosities, accretion rates, and explosion times show,in our models, less variation with the model dimen-sionality than in other studies (see Figure 2). Webelieve this is purely due to the fact that the objec-tive of our study was to obtain models with similarenergetics, rather than to determine the thresholdneutrino luminosity required for the explosion. • We found no evidence for the standing accretionshock instability (SASI) in our models. Our resultsshowed no persistent, oscillatory l = 1 or l = 2mode behavior, considered to be a trademark ofSASI. The absence of SASI in our models is consis-tent with short advective times through the gain re-gion. Instead, the explosions are boosted in multi-dimensions by strong, neutrino- driven convection.At the same time, we would like to stress that SASIcannot be excluded as an important, or even deci-sive, physics mechanism for different parameters ofcollapsing core (e.g., progenitor structure, equationof state, etc.). • We found that the shock front shows a rathermodest degree of asymmetry in our models, r maxs /r mins ≈ .
3. We found tentative evidence be-tween the large-scale structure of convection andthe shock asymmetry. In particular, the shock isless deformed in the case when convection has aricher structure (higher-order modes are present). • We found the neutron star recoil velocities reachingup to ≈
475 km/s, with typical values between 100km/s and 300 km/s. These results are in goodagreement with the previous studies (e.g. Schecket al. 2006). • We introduced and applied new metrics to diagnosethe dynamics and structure of the gain region. Inparticular, we studied the characteristics of the up-flowing material (which we identified as the buoy-ant, convective bubbles) and the solid angle it oc-cupies. Using these measures, we found that thestructure of the gain region on large scales and itsevolution did not significantly differ between 2Dand 3D. However, for the fast contracting mod-els, we found that the mass inside upflows always exceeded 50% of the total gain region, hinting atprompt-like explosions. This situation is qualita-tively different in the case of the slow contractingmodels, where a prolonged period of quasi-steadystate in the gain region can be clearly recognized. Itappears that during that phase the mass in the up-flows fluctuates around the slowly increasing meanvalue with amplitude of several percent, and on atimescale of 10 ms. We speculate that individualdownflows are responsible for the observed fluctu-ations. • Our method of identifying the convective bubblesallowed us to study the statistical properties of thebubble distributions. We found that the solid angleoccupied by bubbles is initially small and progres-sively increases with time (see Figure 15). Thisindicates that the bubbles are small at first andsteadily grow as the evolution progresses. The pro-cess appears to bear some similarity to the bubblemerging process observed in multi-mode Rayleigh-Taylor instabilities. Furthermore, we find a goodcorrelation between the solid angle occupied by ris-ing material and the upflowing mass. This high-lights the convection takes place inside the wholegain region, from the gain radius all the way up tothe shock. • The development and evolution of convection in-side the gain region can also be studied by con-sidering the temporal evolution of the number ofbubbles. We found that their number rapidly in-creased at t ≈
50 ms (see Figure 16), indicatingthe onset of convection. After the convection set in,the number of bubbles progressively decreased withtime. We used the information about the numberof bubbles combined with the solid angle they oc-cupy, and estimated that the convective bubblesare, on average, about 4 times smaller in 3D thanin 2D at the time of the explosion. • We used the decomposition of the total energy fluxto understand the energy transport in the gain re-gion. We found that the work done by expansionand compression in that region differs between 2Dand 3D, and we proposed that this is due to thegreater structural integrity of flow structures in 2D. • The analysis of the components of the total energyflux also provided evidence that the region nearthe gain radius was not in thermal equilibrium atthe beginning of the quasi-steady phase. We foundevidence for strong cooling below and strong heat-ing above the gain radius. Furthermore, we foundthat the buoyant work was positive throughout thegain region in our models, in accordance with theneutrino-driven convection scenario. • We demonstrated that, for the same neutrino lumi-nosities that produce energetic explosions in well-resolved models, the explosion energies dramati-cally decrease once the angular mesh resolution de-creases below 6 ◦ . This implies there exists a mini-mal mesh resolution required for the central engineto efficiently operate. This is expected. This alsohock revival with SN 1987A energetics 23implies that the neutrino-driven convection is thekey ingredient of the explosion mechanism in theenergetic models considered here. We cannot ex-clude the possibility that more physics (e.g. turbu-lence) may emerge at still higher resolutions (i.e.less than about 0 . × r ) to power the explosionmore efficiently than convection. • We found that the the structure of turbulent energyspectra are steeper in our 2D models ( E ∼ l − . )than in our 3D models ( E ∼ l − . ). These resultsshould be taken with caution given the low reso-lution of our models and an estimated Reynoldsnumber on the order of 100. • Analysis of the turbulent Reynolds stresses re-vealed that the post-shock flow is anisotropic inthe radial direction due to buoyant driving. At thesame time, we found the flow in the gain region isisotropic in the lateral directions. • Using a Lagrangian representation of the gain re-gion, we found substantial differences in the en-ergetics of the shocked fluid elements between 2Dand 3D. We identified material near the base of thegain region that experiences stronger heating dur-ing the period when the gain radius shrinks. Wealso found that material accreted during the phasewhen the convection is fully developed undergoesespecially strong heating in 3D. • We found significant differences in the evolution ofresidency times and binding energies of fluid ele-ments inside the gain region between 2D and 3D. In3D, more mass is involved in the process of drivingconvection in the region around the gain radius. Atthe same time, the same amount of positive bind-ing energy is carried by material that entered thegain region during the first 100 ms of evolution asby those which entered the gain region during thefollowing 50 ms. However, in 2D, a smaller portionof the gain region mass participates in driving con-vection, and practically all of the positive bindingenergy is carried by material that entered the gainregion during the first 100 ms.Several possible directions for future research emergefrom the current work. It remains unclear why the con-vective engine is more efficient in three-dimensions. Al-though we were able to quantify the effect, we did notprovide a clear explanation for its origin. We obtainedonly rudimentary information about the properties ofturbulence in our largely underresolved and dominatedby numerical viscosity models. It is clear that to gainmeaningful insight into the role of turbulence in the pro-cess of core-collapse supernova explosion a new genera-tion of far better-resolved models exploiting more effi-cient computational approaches is required. Also, onewould like to evolve those models until the neutron starrecoil velocities saturate. Finally, our current model doesnot include nuclear burning, and therefore we are unableto discuss the nucleosynthetic yields and the composi-tional structure of our models. The above issues will bethe subject of a forthcoming publication. We thank Konstantinos Kifonidis for his comments onthe initial version of the manuscript. The work of TH andTP has been supported by the NSF grant AST-1109113,and TP has been partially supported by the DOE grantDE-FG52-09NA29548. This research used resources ofthe National Energy Research Scientific Computing Cen-ter, which is supported by the Office of Science of the U.S.Department of Energy under Contract No. DE-AC02-05CH11231. Simulations were performed in part usingthe 3Leaf Systems SMP clusters at the Florida StateUniversity High Performance Computing Center, and theDeszno SMP cluster at Jagiellonian University, Cracow,Poland.
REFERENCESAlon, U., Hecht, J., Mukamel, D., & Shvarts, D. 1994, PhysicalReview Letters, 72, 2867Arnett, W. D. 1987, ApJ, 319, 136Blinnikov, S. I. 1999, Astronomy Letters, 25, 359Blondin, J. M., & Mezzacappa, A. 2007, Nature, 445, 58Blondin, J. M., Mezzacappa, A., & DeMarino, C. 2003, ApJ, 584,971Borriello, E., Chakraborty, S., Janka, H.-T., Lisi, E., & Mirizzi,A. 2013, ArXiv e-printsBruenn, S. W. 1993, in Nuclear Physics in the Universe, ed.M. W. Guidry & M. R. Strayer (Bristol: IOP), 31–50Bruenn, S. W., Mezzacappa, A., Hix, W. R., et al. 2013, ApJ,767, L6Brummell, N. H., Clune, T. L., & Toomre, J. 2002, ApJ, 570, 825Burrows, A. 2013, Reviews of Modern Physics, 85, 245Burrows, A., & Goshy, J. 1993, ApJ, 416, L75Chugai, N. N. 1988, Pis’ma Astronomicheskii Zhurnal, 14, 1079Colella, P., & Glaz, H. M. 1985, J. Comput. Phys., 59, 264Colella, P., & Woodward, P. R. 1984, J. Comput. Phys., 54, 174Couch, S. M. 2013, ApJ, 775, 35Dolence, J. C., Burrows, A., Murphy, J. W., & Nordhaus, J.2013, ApJ, 765, 110Einfeldt, B. 1988, SIAM J. Num. Anal., 25, 294Endeve, E., Cardall, C. Y., Budiardja, R. D., et al. 2012, ApJ,751, 26Foglizzo, T. 2009, ApJ, 694, 820Foglizzo, T., Scheck, L., & Janka, H.-T. 2006, ApJ, 652, 1436Frisch, U. 1995, Turbulence: the legacy of A.N. KolmogorovFryer, C. L., & Warren, M. S. 2002, ApJ, 574, L65Fryer, C. L., & Young, P. A. 2007, ApJ, 659, 1438Gawryszczak, A., Guzman, J., Plewa, T., & Kifonidis, K. 2010,A&A, 521, A38Hammer, N. J., Janka, H.-T., & M¨uller, E. 2010, ApJ, 714, 1371Hanke, F., Marek, A., M¨uller, B., & Janka, H.-T. 2012, ApJ, 755,138Herant, M., Benz, W., Hix, W. R., Fryer, C. L., & Colgate, S. A.1994, ApJ, 435, 339Hurlburt, N. E., Toomre, J., & Massaguer, J. M. 1986, ApJ, 311,563Imshennik, V. S., & Popov, D. V. 1992, AZh, 69, 497Janka, H.-T., Buras, R., Kifonidis, K., Plewa, T., & Rampp, M.2003, in From Twilight to Highlight: The Physics ofSupernovae, ed. W. Hillebrandt & B. Leibundgut, 39Janka, H.-T., Hanke, F., H¨udepohl, L., et al. 2012, Progress ofTheoretical and Experimental Physics, 2012, 010000Janka, H.-T., & Mueller, E. 1996, A&A, 306, 167Kifonidis, K., Plewa, T., Janka, H.-T., & M¨uller, E. 2003, A&A,408, 621Kifonidis, K., Plewa, T., Scheck, L., Janka, H.-T., & M¨uller, E.2006, A&A, 453, 661Laming, J. M. 2007, ApJ, 659, 1449Lattimer, J. M., & Swesty, D. F. 1991, Nuclear Physics A, 535,331Marek, A., Dimmelmeier, H., Janka, H.-T., M¨uller, E., & Buras,R. 2006, A&A, 445, 273Mezzacappa, A., Calder, A. C., Bruenn, S. W., et al. 1998, ApJ,495, 911
Miles, A. R. 2004, Physics of Plasmas, 11, 5140Moc´ak, M., M¨uller, E., Weiss, A., & Kifonidis, K. 2009, A&A,501, 659Monin, A. S., & I’Aglom, A. M. 1971, Statistical fluid mechanics;mechanics of turbulenceM¨uller, B., Janka, H.-T., & Heger, A. 2012, ApJ, 761, 72Murphy, J. W., & Burrows, A. 2008, ApJ, 688, 1159Murphy, J. W., & Meakin, C. 2011, ApJ, 742, 74Nordhaus, J., Burrows, A., Almgren, A., & Bell, J. 2010, ApJ,720, 694Ohnishi, N., Kotake, K., & Yamada, S. 2006, ApJ, 641, 1018Ott, C. D., Abdikamalov, E., M¨osta, P., et al. 2013, ApJ, 768, 115Pumo, M. L., & Zampieri, L. 2011, ApJ, 741, 41Rampp, M., & Janka, H.-T. 2002, A&A, 396, 361Scheck, L. 2007, Phd thesis, Technische Universit¨at M¨unchenScheck, L., Janka, H.-T., Foglizzo, T., & Kifonidis, K. 2008,A&A, 477, 931 Scheck, L., Kifonidis, K., Janka, H.-T., & M¨uller, E. 2006, A&A,457, 963Shigeyama, T., & Nomoto, K. 1990, ApJ, 360, 242Strang, G. 1968, SIAM J. Numer. Anal., 5, 506Takiwaki, T., Kotake, K., & Suwa, Y. 2013, ArXiv e-printsUtrobin, V. 1993, A&A, 270, 249Utrobin, V. P. 2004, Astronomy Letters, 30, 293—. 2005, Astronomy Letters, 31, 806Utrobin, V. P., & Chugai, N. N. 2005, A&A, 441, 271Wongwathanarat, A., Janka, H.-T., & M¨uller, E. 2013, A&A, 552,A126Woosley, S. E. 1988, ApJ, 330, 218 hock revival with SN 1987A energetics 25
Figure 23.
Analysis results for models M194J and M187D. (left panel) Estimated numerical kinematic viscosity. (right panel) EstimatedReynolds number as a function of particle separation distance. We may expect our models to reflect flows with Reynolds numbers on theorder of 10s or 100s. APPENDIX
ESTIMATING NUMERICAL VISCOSITY AND REYNOLDS NUMBER
In order to estimate the numerical viscosity and the numerical Reynolds number, we use the information providedby the Lagrangian tracer particles. Consequently, in our approach the velocity field is represented at discrete pointsby tracer particles. We define the longitudinal Eulerian velocity structure function of order- p as S Lp ( r sep ) = (cid:104) [( u ( x ) − u ( x + r sep ˆ r )) · ˆ r ] p (cid:105) , (A1)where x is the particle position, u is the particle velocity, r sep is the separation between particles, and ˆ r is the particledisplacement unit vector. The averaging operator, (cid:104)·(cid:105) , is applied over a given r sep at a specific model time.Kolmogorov’s turbulence theory (Monin & I’Aglom 1971; Frisch 1995), assuming negligible intermittence, states thatfor p = 2, the structure function, S L , is related to the average specific turbulent dissipation, (cid:104) (cid:15) (cid:105) , via S L ( r sep ) = C (cid:104) (cid:15) (cid:105) / r / sep , (A2)where C is a constant on the order of 1. This result allows one to estimate the turbulent dissipation (physical +numerical) present in a numerical simulation. Within the same framework, the local, specific turbulent dissipation isrelated to the average strain rate as (cid:15) = 2 ν (cid:104) S : S (cid:105) , (A3)where ν is the kinematic viscosity and S is the strain rate tensor.In order to estimate the numerical Reynolds number in our models, we assume that the local turbulent dissipationis representative of the average dissipation rate, (cid:15) ≈ (cid:104) (cid:15) (cid:105) . This allows us to compute the kinematic viscosity as ν = S L / C / (cid:104) S : S (cid:105) r sep . (A4)Using r sep as our characteristic length scale and approximating the characteristic velocity in the field by (cid:112) S L , weestimate the Reynolds number as a function of separation distance via Re = r sep (cid:112) S L ν . (A5)Figure 23 shows the numerical kinematic viscosity and numerical Reynolds number estimates for the 2D SC modelM194J and the 3D SC model M187D at t = 150 ms as a function of particle separation distance. We find numericalviscosities on the order of 1 × cm s − for the 2D model and 1 × cm s − for the 3D model. At this timein the model evolution, the average shock radius is about 500 km, providing an upper limit on structure size. Onthat scale, the numerical Reynolds number is ≈
200 in 2D and ≈≈