Massive Computation for Understanding Core-Collapse Supernova Explosions
CC. D. OTT, MANUSCRIPT FOR COMPUTING IN SCIENCE AND ENGINEERING (CISE), DATED: SEPTEMBER 4, 2018 1
Massive Computation for UnderstandingCore-Collapse Supernova Explosions
Christian D. Ott
TAPIR, Walter Burke Institute for Theoretical Physics, California Institute of Technology, Pasadena, CA, 91125, USA Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto, Japan
How do massive stars explode? Progress toward the answer is driven by increases in compute power. Petascale supercomputersare enabling detailed three-dimensional simulations of core-collapse supernovae. These are elucidating the role of fluid instabilities,turbulence, and magnetic field amplification in supernova engines.
Index Terms —Supernovae, neutron stars, gravitational collapse
I. I
NTRODUCTION
Core-collapse supernova explosions come from stars moremassive than ∼ − times the mass of the Sun. Tencore-collapse supernovae explode per second in the universe,automated astronomical surveys discover multiple per night,and one or two explode per century in the Milky Way. Core-collapse supernovae outshine entire galaxies in photons forweeks and output more power in neutrinos than the combinedlight output of all other stars in the universe, for tens ofseconds. These explosions pollute the interstellar mediumwith the ashes of thermonuclear fusion. From these elements,planets form and life is made. Supernova shock waves stir theinterstellar gas, trigger or shut off the formation of new stars,and eject hot gas from galaxies. At their centers, a stronglygravitating compact remnant, a neutron star or a black hole,is formed.As the name alludes, the explosion is preceded by collapseof a stellar core. At the end of its life, a massive star hasa core composed mostly of iron-group nuclei. The core issurrounded by an onion-skin structure of shells dominated bysuccessively lighter elements. Nuclear fusion is still ongoingin the shells, but the iron core is inert. The electrons in the coreare relativistic and degenerate. They provide the lion’s share ofthe pressure support stabilizing the core against gravitationalcollapse. In this, the iron core is very similar to a whitedwarf star, the end product of low-mass stellar evolution.Once the iron core exceeds its maximum mass (the so-calledeffective Chandrasekhar mass of ∼ . − solar masses [ M (cid:12) ]),gravitational instability sets in. Within a few tenths of asecond, the inner core collapses from a central density of ∼ g cm − to a density comparable to that in an atomicnucleus ( (cid:38) . × g cm − ). There, the repulsive part ofthe nuclear force causes a stiffening of the equation of state(EOS; the pressure–density relationship). The inner core firstovershoots nuclear density, then rebounds (“bounces”) intothe still collapsing outer core. The inner core then stabilizesand forms the inner regions of the newborn protoneutronstar. The hydrodynamic supernova shock is created at theinterface of inner and outer core. First, the shock movesoutward dynamically. It then quickly loses energy by work H He Iron CoreRed Supergiant (not drawn to scale)
SiC,O
R≈10 km InnerCore1.5−2 M ⊙ ≈ k m Accretion ≈ k m Core Collapse toProtoneutron Star(PNS)Stalled Shock
PNS ν ν ν ννννννν M • M • M • M • ShockrevivedShock notrevived © A n g l o - A u s t r a l i a n O b s e r v a t o r y Core-CollapseSupernova ExplosionBlack HoleFormationτ ≈ 1− fewSeconds
Fig. 1. Schematic of core collapse and its simplest outcomes. The imageshows SN 1987A, which exploded in the Large Magellanic Cloud. a r X i v : . [ a s t r o - ph . H E ] A ug . D. OTT, MANUSCRIPT FOR COMPUTING IN SCIENCE AND ENGINEERING (CISE), DATED: SEPTEMBER 4, 2018 2 done breaking up infalling iron-group nuclei into neutrons,protons, and alpha particles. The copious emission of neutrinosfrom the hot ( T ∼
10 MeV (cid:39) K ) gas further reducesenergy and pressure behind the shock. The shock stalls andturns into an accretion shock: the ram pressure of accretion ofthe star’s outer core balances the pressure behind the shock.The supernova mechanism must revive the stalled shockto drive a successful core-collapse supernova explosion. De-pending on the structure of the progenitor star, this mustoccur within one to a few seconds of core bounce. Otherwise,continuing accretion pushes the protoneutron star over itsmaximum mass ( ∼ − M (cid:12) ), which results in the formationof a black hole and no supernova explosion.If the shock is successfully revived, it must travel throughthe outer core and the stellar envelope before it breaks out ofthe star and creates the spectacular explosive display observedby astronomers on Earth. This may take more than a day fora red supergiant star (e.g., like Betelgeuse, a ∼ M (cid:12) star inthe constellation Orion) or just tens of seconds for a star thathas been stripped of its extended hydrogen-rich envelope bya strong stellar wind or mass exchange with a companion starin a binary system.The photons observed by astronomers are emitted extremelyfar from the central regions. They carry information on theoverall energetics, the explosion geometry, and on the productsof explosive nuclear burning that is triggered by the passingshock wave. They can, however, only provide weak constraintson the inner workings of the supernova. Direct observationalinformation on the supernova mechanism can be gained onlyfrom neutrinos and gravitational waves that are emitted di-rectly in the supernova core. Detailed computational modelsare required for gaining theoretical insight and for makingpredictions that can be contrasted with future neutrino andgravitational-wave observations from the next core-collapsesupernova in the Milky Way.II. S UPERNOVA E NERGETICS AND M ECHANISMS
Core-collapse supernovae are “gravity bombs.” The energyreservoir from which any explosion mechanism must drawis the gravitational energy released in the collapse of theiron core to a neutron star: ∼ × erg ( × J ), amass-energy equivalent of ∼ . M (cid:12) c . A fraction of thistremendous energy is stored initially as heat (and rotationalkinetic energy) in the protoneutron star and the rest comesfrom its subsequent contraction. Astronomical observations,on the other hand, show the typical core-collapse supernovaexplosion energy to be in the range − erg . Hypernova explosions may have up to erg , but they make up (cid:46) (cid:38) ( (cid:38) in the hypernova case) of the available energy over O (10) s as the protoneutron star cools and contracts. This wasfirst theorized and then later observationally confirmed withthe detection of neutrinos from SN 1987A, the most recentcore-collapse supernova in the Milky Way vicinity. Fig. 2. Volume rendering of the specific entropy in the core of a neutrino-driven core-collapse supernova at the onset of explosion. Based on the 3Dgeneral-relativistic simulations of [1] and rendered by Steve Drasco (Cal PolySan Luis Obispo). Specific entropy is a preferred quantity for visualization,since in the core of a supernova, it typically ranges from ∼ ∼
20 unitsof Boltzmann’s constant k B per baryon. Shown is the large-scale asymmetricshock front and a layer of hot expanding plumes behind it. The physical scaleis roughly ×
400 km . Since neutrinos dominate the energy transport through thesupernova, they might quite naturally have something to dowith the explosion mechanism. The neutrino mechanism , inits current form, was proposed by Bethe & Wilson [2]. Inthis mechanism, a fraction ( ∼ ) of the outgoing electronneutrinos and antineutrinos is absorbed in a layer betweenprotoneutron star and the stalled shock. In the simplest picture,this neutrino heating increases the thermal pressure behind thestalled shock. Consequently, the dynamical pressure balanceat the accretion shock is violated and a runaway explosion islaunched.The neutrino mechanism fails in spherical symmetry (1D,e.g., [3]), but is very promising in multiple dimensions (ax-isymmetry [2D], 3D). This is due largely to multi-D hydrody-namic instabilities that break spherical symmetry (see Figure 2for an example), increase the neutrino mechanism’s efficiency,and facilitate explosion. I will discuss this in more detail laterin this article. The neutrino mechanism is presently favored asthe mechanism driving most core-collapse supernova explo-sions (see [3] for a recent review).Despite its overall promise, the neutrino mechanism is veryinefficient. Only (cid:46) of the outgoing total electron neutrinoand antineutrino luminosity is deposited behind the stalledshock at any moment and much of this deposition is lost again . D. OTT, MANUSCRIPT FOR COMPUTING IN SCIENCE AND ENGINEERING (CISE), DATED: SEPTEMBER 4, 2018 3 Fig. 3. Volume rendering of the specific entropy in the core of a magne-torotational core-collapse supernova. Bluish colors indicate low entropy, redcolors high entropy, and green and yellow intermediate entropy. The verticalis the axis of rotation and shown is a region of ∼ ×
800 km . The ultra-strong toroidal magnetic field surrounding the the protoneutron star pushes hotplasma out along the rotation axis. The distorted, double-lobe structure is dueto an MHD kink instability akin those seen in Tokamak fusion experiments.Used with permission from M¨osta et al. as heated gas flows down, leaves the heating region, and settlesonto the protoneutron. The neutrino mechanism may (barely)be able to power ordinary core-collapse supernovae, but itcannot deliver hypernova explosion energies or account forgamma-ray bursts.An alternative mechanism that may be part of the explana-tion for such extreme events is the magnetorotational mecha-nism , first suggested by Bisnovatyi-Kogan [5] and LeBlanc &Wilson [6]. In its modern form, a very rapidly spinning corecollapses to a protoneutron star with a spin period of only ∼ . Its core is expected to be spinning uniformly,but its outer regions will be extremely differentially rotating.These are ideal conditions for the magnetorotational instability (MRI, [7]) to operate, amplify any seed magnetic field, anddrive magnetohydrodynamic (MHD) turbulence. If a dynamoprocess is present, an ultra-strong large-scale (globally or-dered) magnetic field is built up. This makes the protoneutronstar a protomagnetar . Provided this occurs, magnetic pressuregradients and hoop stresses could lead to outflows along theaxis of rotation. The MRI’s fastest growing mode has a smallwavelength and is extremely difficult to resolve numerically. Because of this, all simulations of the magnetorotationalmechanism to date have simply made the assumption that acombination of MRI and dynamo is operating. They then ad-hoc imposed a strong large-scale field as an initial condition.In 2D simulations, collimated jets develop along the axis ofrotation. In 3D, the jets are unstable and a more complicatedexplosion geometry develops [4], as shown in Figure 3. Nev-ertheless, even in 3D, an energetic explosion could potentiallybe powered.The magnetorotational mechanism requires one specialproperty of the progenitor star: rapid core rotation. Presently,stellar evolution theory suggests that the cores of most massivestars should be slowly spinning. However, there may beexceptions of rapidly spinning cores at just about the rightoccurrence rate to explain hypernovae and long gamma-raybursts.Besides the neutrino mechanism and the magnetorotationalmechanism, a number of other explosion mechanisms havebeen proposed. I direct the interested reader to the moreextensive review by [3].III. A M ULTI -S CALE , M
ULTI -P HYSICS ,M ULTI -D IMENSIONAL C OMPUTATIONAL C HALLENGE
The core-collapse supernova problem is highly complex,inherently non-linear, and involves many branches of (as-tro)physics. Only limited progress can be made with analyticor perturbative methods. Computational simulation is a pow-erful means for gaining theoretical insight and for making pre-dictions that could be tested with astronomical observations ofneutrinos, gravitational waves, and electromagnetic radiation.Core-collapse supernova simulations are time evolution sim-ulations – starting from initial conditions, the matter, radiation,and gravitational fields are evolved in time. In the case oftime-explicit evolution, the numerical timestep is limited bycausality, controlled by the speed of sound in Newtoniansimulations, and the speed of light in general-relativistic simu-lations. Because of this, an increase in the spatial resolution bya factor of two corresponds to a decrease in the time step bya factor of two. Hence, in a 3D simulation, the computationalcost scales with the fourth power of resolution.
A. Multi Scale
Taking the red supergiant in Figure 1 as an example, acomplete core-collapse supernova simulation that follows theshock to the stellar surface, would have to cover dynamics on aphysical scale from ∼ km (stellar radius) down to ∼ . (the typical scale over which structure and thermodynamics ofthe protoneutron star change). These ten orders of magnitudein spatial scale are daunting. In practice, reviving the shockand tracking its propagation to the surface can be treated as(almost) independent problems. If our interest is on the shockrevival mechanism, we need to include the inner ∼ ,
000 km of the star. Since information about core collapse is commu-nicated to overlying layers with the speed of sound, stellarmaterial at greater radii will not “know” that core collapsehas occurred before it is hit by the revived expanding shock. . D. OTT, MANUSCRIPT FOR COMPUTING IN SCIENCE AND ENGINEERING (CISE), DATED: SEPTEMBER 4, 2018 4
Even with only five decades in spatial scale, some formof grid refinement or adaptivity is called for: a 3D finite-difference grid with an extent of ,
000 km symmetric aboutthe origin with uniform . cell size would require 57PB of RAM to store a single double precision variable.Many tens to hundreds of 3D variables are required. Suchhigh uniform resolution is not only currently impossible butalso unnecessary. Most of the resolution is needed near theprotoneutron star and in the region behind the stalled shock.The near-free-fall collapse of the outer core can be simulatedwith much lower resolution.Because of the broad range of physics involved (see below)and the limited available compute power, early core-collapsesupernova simulations were spherically symmetric (1D). 1Dsimulations often employ a Lagrangian, comoving mass co-ordinate discretization. This grid can be set up to providejust the right resolution where and when needed or can bedynamically re-zoned (an adaptive mesh refinement [AMR]technique). Other 1D codes discretize in the Eulerian frameand use a fixed grid whose cells are radially stretched usinggeometric progression.In 2D simulations, Eulerian, geometrically-spaced fixedspherical grids are the norm, but some codes use cylindricalcoordinates and AMR. Spherical grids, already in 2D, sufferfrom a coordinate singularity at the axis that can lead tonumerical artifacts. In 3D, they become even more difficult tohandle and their focusing grid lines impose a severe timestepconstraint near the origin. Some 3D codes still use a spher-ical grid, while many others employ Cartesian AMR grids.Recent innovative approaches use so-called multi-block gridswith multiple curvilinear touching or overlapping logicallyCartesian “cubed-sphere” grids (e.g., [8]). B. Multi Physics
Core-collapse supernovae are very rich in physics. Allfundamental forces are involved and essential to the corecollapse phenomenon. These forces are probed under condi-tions that are impossible (or exceedingly difficult) to create inearthbound laboratories.Gravity drives the collapse and provides the energy reser-voir. It is so strong near the protoneutron star that general rel-ativity becomes important and its Newtonian description doesnot suffice. The electromagnetic force describes the interactionof the dense, hot magnetized, perfectly conducting plasmaand the photons that provide thermal pressure and make thesupernova light. The weak force governs the interactions ofneutrinos and the strong (nuclear) force is essential in thenuclear EOS and nuclear reactions.All this physics occurs at the microscopic, per particlelevel. Fortunately, the continuum assumption holds, allowingus to describe core-collapse supernovae on a macroscopic scaleby a coupled set of systems of non-linear partial differentialequations (PDEs): • (Magneto)hydrodynamics (MHD). The stellar plasma isin local thermodynamic equilibrium, essentially perfectlyconducting, and essentially inviscid (though neutrinosmay provide some shear viscosity in the protoneutron
30 km60 km120 km150 km240 km
Fig. 4. Map projections of the momentum-space neutrino radiation field(for ν e at an energy of 16.3 MeV) going outward radially (from top tobottom) on the equator of a supernova core. Generated using the simulationresults of [9]. Inside the protoneutron star ( R (cid:46)
30 km ) neutrinos andmatter are in equilibrium and the radiation field is isotropic. It becomesmore and more forward peaked as the neutrinos decouple and become freestreaming. Handling the transition from slow diffusion to free streamingcorrectly requires angle-dependent radiation transport, which is a 6+1 Dproblem and computationally extremely challenging. star). The ideal, inviscid MHD approximation is appro-priate under these conditions. The MHD equations arehyperbolic and can be written in flux-conservative formwith source terms that do not include derivatives of theMHD variables. They are typically solved with standardtime-explicit high-resolution shock capturing methodsthat exploit the characteristic structure of the equations(e.g., [10]). Special attention must be paid to preservingthe divergence-free property of the magnetic field. TheMHD equations require an EOS as a closure (see below).Unless ultra-strong ( B (cid:38) G ), magnetic fieldshave little effect on the supernova dynamics and thus arefrequently neglected. Since strong gravity and velocitiesup to a few tenths of the speed of light are involved, theMHD equations are best solved in a general-relativisticformulation. General-relativistic MHD is computationallyparticularly expensive, because the conserved variablesare not the primitive variables (density, internal energy /temperature, velocity, chemical composition). The latterare needed for the EOS and enter flux terms. After . D. OTT, MANUSCRIPT FOR COMPUTING IN SCIENCE AND ENGINEERING (CISE), DATED: SEPTEMBER 4, 2018 5 each update, they must be recovered from the conservedvariables via multi-dimensional root finding. • Gravity . Deviations in the strength of the gravitationalacceleration between Newtonian and general-relativisticgravity are small in the precollapse core, but become oforder − in the protoneutron star phase. In the caseof black hole formation, Newtonian physics breaks downcompletely. General relativistic gravity is included atvarying levels in simulations. Some neglect it completelyand solve the linear elliptic Newtonian Poisson equationto compute the gravitational potential. This is done usingdirect multigrid methods or integral multipole expansionmethods. Some codes modify the monopole term in thelatter approach to approximate general relativistic effects.Including full general relativity is more challenging,in particular in 2D and 3D, since there general relativityhas radiative degrees of freedom (gravitational waves).An entire subfield of gravitational physics, numericalrelativity , spent nearly five decades looking for waysto solve Einstein’s equations on computers (see [11]for a comprehensive introduction). In general relativity,changes in the gravitational field propagate at the speedof light. Hence, time evolution equations must be solved.This is done by splitting 4D spacetime into 3D spatialslices that are evolved in the time direction. In thesimplest way of writing the equations (the so-calledArnowitt-Deser-Misner [ADM] formulation), they forma system of 12 partial differential evolution equations,4 gauge variables that must be specified (and evolvedin time or recalculated on each slice), and 4 elliptic constraint equations without time derivatives. The ADMformulation has poor numerical stability properties. Theselead to violations of the constraint equations and numeri-cal instabilities that make long-term evolution impossible.It took until the 2000s for numerical relativity to findformulations of Einstein’s equations and gauge choicesthat together lead to stable long-term evolutions. Insome cases, well-posedness and strong or symmetrichyperbolicity can be proven. The equations are typicallyevolved time-explicitly with straightforward high-order(fourth and higher) finite difference schemes or withmulti-domain pseudospectral methods.Since numerical relativity only recently became appli-cable to astrophysical simulations, very few core-collapsesupernova codes are fully general relativistic at thispoint [1], [12]. The fully general-relativistic approach ismuch more memory and FLOP intensive than solvingthe Newtonian Poisson equation. Its advantage in large-scale computations, however, is the hyperbolic natureof the equations, which does not require global matrixinversions or summations and thus is advantageous forthe parallel scaling of the algorithm. • Neutrino Transport and Neutrino-Matter Interac-tions.
Neutrinos move at the speed of light (the very smallneutrino masses are neglected) and can travel macro-scopic distances between interactions. Therefore, theymust be treated as non-equilibrium radiation. Radiation transport is closely related to kinetic theory’s Boltzmannequation. It describes the phase-space evolution of theneutrino distribution function or, in radiation transportterminology, their specific intensity. This is a 6+1 Dproblem: 3 spatial dimensions, neutrino energy, and twomomentum space propagation angles in addition to time.The angles describe the directions from which neutrinosare coming and where they are going at a given spatialcoordinate. In addition, the transport equation must besolved separately for multiple neutrino species: electronneutrinos, electron antineutrinos, and heavy-lepton ( µ , τ )neutrinos and antineutrinos.Figure 4 shows map projections of the momentumspace angular neutrino distribution at different radii ina supernova core. In the dense protoneutron star, neutri-nos are trapped and in equilibrium with matter. Theirradiation field is isotropic. They gradually diffuse outand decouple from matter at the neutrinosphere (theneutrino equivalent of the photosphere). This decouplingis gradual and marked by the transition of the angulardistribution into the forward (radial) direction. In theouter decoupling region, neutrino heating is expected tooccur and the heating rates are sensitive to the angulardistribution of the radiation field (cf. [9]). Eventually, atradii of a few hundred kilometers, the neutrinos have fullydecoupled and are free streaming. Neutrino interactionswith matter (and thus the decoupling process) are verysensitive to neutrino energy, since weak-interaction cross-sections scale with the square of the neutrino energy.This is why neutrino transport needs to be multi-group ,with typically a minimum of − energy groupscovering supernova neutrino energies of −O (100) MeV .Typical mean energies of electron neutrinos are around −
30 MeV . Energy exchanges between matter andradiation occur via the collision terms in the Boltz-mann equation. These are stiff sources/sinks that mustbe handled time-implicitly with (local) backward-Eulermethods. The neutrino energy bins are coupled through(1) frame-dependent energy shifts since the material neu-trinos interact with is moving, (2) gravitational redshift,and (3) energy transfer in scatterings off of electrons andnucleons. Neutrino-matter interaction rates are usuallyprecomputed and stored in dense multi-D tables withinwhich simulations interpolate.Full 6+1 D general-relativistic Boltzmann neutrino-radiation hydrodynamics is exceedingly challenging andhas so far not been possible to included in core-collapsesupernova simulations. 3+1 D (1D in space, 2D inmomentum space) (e.g., [13]), 5+1 D (2D in space,3D in momentum space) simulations [9] and static 6Dsimulations [14] have been carried out.Most (spatially) multi-D simulations treat neutrinotransport in some dimensionally-reduced approximation.The most common is an expansion of the radiation fieldinto angular moments. The n -th moment of this expansionrequires information about the ( n + 1) -th moment (andin some cases also about the ( n + 2) -th moment). Thisnecessitates a closure relation for the moment at which . D. OTT, MANUSCRIPT FOR COMPUTING IN SCIENCE AND ENGINEERING (CISE), DATED: SEPTEMBER 4, 2018 6 the expansion is truncated. Multi-group flux-limited dif-fusion evolves the -th moment (the radiation energydensity). The flux limiter is the closure that interpolatesbetween diffusion and free streaming. The disadvantagesof this method are its very diffusive nature that washesout spatial variations of the radiation field, its sensitivityto the choice of flux limiter, and the need for time-implicitintegration (involving global matrix inversion) due to thestability properties of the parabolic diffusion equation.Two-moment transport is the next better approximation.It solves equations for the radiation energy density andmomentum (i.e. the radiative flux) and requires a closurethat describes the radiation pressure tensor (also knownas the Eddington tensor). This closure can be analytic andbased on the local values of energy density and flux (theM1 approximation). Alternatively, some codes compute aglobal closure based on the solution of a simplified, time-independent Boltzmann equation. The major advantageof the two-moment approximation is that its advectionterms are hyperbolic and can be handled with standardtime-explicit finite-volume methods of computational hy-drodynamics and only the local collision terms need time-implicit updates.There are now implementations of multi-group two-moment neutrino radiation-hydrodynamics in multiple2D/3D core-collapse supernova simulation codes (e.g.,[12], [15], [16]). This method may be sufficiently closeto the full Boltzmann solution (in particular if a globalclosure is used) and appears to be the way towardmassively-parallel long-term 3D core-collapse supernovasimulations. • Neutrino Oscillations.
Neutrinos have mass and can os-cillate between flavors. The oscillations occur in vacuum,but can also be mediated by neutrino-electron scattering(the Mikheyev-Smirnov-Wolfenstein [MSW] effect) andneutrino-neutrino scattering. Neutrino oscillations dependon neutrino mixing parameters and on the neutrino masseigenstates (the magnitudes of the mass differences areknown, but not their signs). Observation of neutrinosfrom the next galactic core-collapse supernova could helpconstrain the neutrino mass hierarchy (see the recentreview by [17]).MSW oscillations occur in the stellar envelope. Theyare important for the neutrino signal observed in detectorson Earth, but they cannot influence the explosion itself.The self-induced (via neutrino-neutrino scattering) oscil-lations, however, occur at the extreme neutrino densitiesnear the core. They offer a rich phenomenology thatincludes collective oscillation behavior of neutrinos (seethe review in [17]). The jury is still out on their potentialinfluence on the explosion mechanism.Collective neutrino oscillation calculations (essentiallysolving coupled Schr¨odinger-like equations) are compu-tationally intensive [17]. They are currently performedindependently of core-collapse supernova simulations anddo not take into account feedback on the stellar plasma.Fully understanding collective oscillations and their im- pact on the supernova mechanism will quite likely requirethat neutrino oscillations, transport, and neutrino-matterinteractions are solved for together in a quantum-kineticapproach [18]. • Equation of State and Nuclear Reactions.
The EOSis essential for the (M)HD part of the problem andfor updating the matter thermodynamics after neutrino-matter interactions. Baryons (proton, neutrons, alpha par-ticles, heavy nuclei), electrons, positrons, and photonscontribute to the EOS. Neutrino momentum transfer con-tributes an effective pressure that is taken into accountseparately since neutrinos are not everywhere in localthermodynamic equilibrium with the stellar plasma. Indifferent parts of the star, different EOS physics applies.At low densities and temperatures below ∼ . ( ∼ × K ), nuclear reactions are too slow to reachnuclear statistical equilibrium. In this regime, the massfractions of the various heavy nuclei (isotopes, in thefollowing) must be tracked explicitly. As the core col-lapses, the gas heats up and nuclear burning must betracked with a nuclear reaction network, a stiff system ofODEs. Solving the reaction network requires the inver-sion of sparse matrices at each grid point. Depending onthe number of isotopes tracked (ranging, typically from O (10) to O (100) ), nuclear burning can be a significantcontributor to the overall computational cost of a simu-lation. The EOS in the burning regime is simple, sinceall isotopes can essentially be treated as non-interactingideal Boltzmann gases. Often, corrections for Coulombinteractions are included. Photons and electrons/positronscan be treated everywhere as ideal Bose and Fermi gases,respectively. Since electrons will be partially or com-pletely degenerate, computing the electron/positron EOSinvolves the FLOP-intensive solution of Fermi integrals.Because of this, their EOS is often included in tabulatedform.At temperatures above ∼ . , nuclear statisticalequilibrium holds. This greatly simplifies things, sincenow the electron fraction Y e (number of electrons perbaryon; because of macroscopic charge neutrality, Y e is equal to Y p , the number fraction of protons) is theonly compositional variable. The mass fractions of allother baryonic species can be obtained by solving Saha-like equations for compositional equilibrium. At densitiesbelow ∼ − g cm − the baryons can still betreated as ideal Boltzmann gases (but including Coulombcorrections).The nuclear force becomes relevant at densities nearand above − g cm − . It is an effective quantummany-body interaction of the strong force and its detailedproperties are presently not known. Under supernovaconditions, matter will be in NSE in the nuclear regimeand the EOS is a function of density, temperature, and Y e . Starting from a nuclear force model, an EOS can beobtained in multiple ways (see the [19] for an overviewdiscussion), including direct Hartree-Fock many-bodycalculations, mean field models, or phenomenological . D. OTT, MANUSCRIPT FOR COMPUTING IN SCIENCE AND ENGINEERING (CISE), DATED: SEPTEMBER 4, 2018 7 Hydrodynamics / MHDNeutrino Transport and InteractionsGravityEquation of State / Nuclear ReactionsSimulationFramework
AMRMemoryCouplingSchedulingCommunicationI/O
Core-Collapse Supernova Simulation Components
Fig. 5. Multi-physics modules of core-collapse supernova simulation codes.The simulation framework provides parallelization, I/O, execution scheduling,AMR, and memory management. models (e.g., the liquid-drop model). Typically, the mini-mum of the Helmholtz free energy is sought and all ther-modynamic variables are obtained from derivatives of thefree energy. In most cases, EOS calculations are too timeconsuming to be performed during a simulation. As in thecase of the electron/positron EOS, large ( (cid:38)
200 MB; mustbe stored by each
MPI process), densely spaced nuclearEOS tables are precomputed and simulations efficientlyinterpolate in (log ρ, log T, Y e ) to obtain thermodynamicand compositional information. C. Effects of Multidimensionality
Stars are, at zeroth order, gas spheres. It is thus naturalto start with assuming spherical symmetry in simulations –in particular given the very limited compute power availableto the pioneers of supernova simulations. After decades ofwork, it appears now clear that detailed spherically symmetricsimulations robustly fail at producing explosions for stars thatare observed to explode in nature. Spherical symmetry itselfmay be the culprit, since symmetry is clearly broken in core-collapse supernovae:( i ) Observations show that neutron stars receive “birthkicks” giving them typical velocities of O (100) km s − withrespect to the center of mass of their progenitors. The mostlikely and straightforward explanation for these kicks arehighly asymmetric explosions leading to neutron star recoilowing to momentum conservation.( ii ) Deep observations of supernova remnants show thatthe innermost supernova ejecta exhibit low-mode asphericitysimilar to the geometry of the shock front shown in Figure 2.( iii ) Analytic considerations and also 1D core-collapsesimulations show that the protoneutron star and the regionbehind the stalled shock where neutrino heating takes placeare both unstable to buoyant convection, which always leadsto the breaking of spherical symmetry.( iv ) Rotation and magnetic fields naturally break sphericalsymmetry. Observations of young pulsars show that someneutron stars must be born with rotation periods of order
10 milliseconds . Magnetars may be born with even shorterspin periods if their magnetic field is derived from rapiddifferential rotation. ( v ) Multi-D simulations of the violent nuclear burning inthe shells overlying the iron core show that large-scale devia-tions from sphericity develop that couple into the precollapseiron core via the excitation of non-radial pulsations [20]. Thesecreate perturbations from which convection will grow aftercore bounce.Given the above, multi-D simulations are essential forstudying the dynamics of the supernova engine.The rapid increase of compute power since the early1990s has facilitated increasingly detailed 2D radiation-hydrodynamics simulations over the past two and a halfdecades. 3D simulations with simplified neutrino treatmentshave been carried out since the early 2000s. The first 3Dneutrino radiation-hydrodynamics simulations have becomepossible only in the past few years, thanks to the computepower of large petascale systems like US NSF/NCSA BlueWaters, US DOE/ORNL Titan or the Japanese K computer.IV. C ORE -C OLLAPSE S UPERNOVA S IMULATION C ODES
Many 1D codes exist, some are no longer in use, and oneis open source and free to download (http://GR1Dcode.org).There are ∼ (depending on how one counts) multi-D core-collapse supernova simulation codes in the community. Many,in particular the 3D codes, follow the design encapsulatedby Figure 5. They employ a simulation framework (e.g., FLASH , http://flash.uchicago.edu/site/flashcode/ or
Cactus http://cactuscode.org) that handles domain decomposition,message passing, memory management, AMR, coupling ofdifferent physics components, execution scheduling, and I/O.Given the tremendous memory requirement and FLOP-consumption of the core-collapse supernova problem, thesecodes are massively parallel and employ both node-local
OpenMP and inter-node
MPI parallelization. All current codesfollow a data-parallel paradigm with monolithic sequentialscheduling. This limits scaling, can create load imbalanceswith AMR, and makes the use of GPU/MIC accelerators chal-lenging, since communication latencies between acceleratorand CPU block execution in the current paradigm.The Caltech
Zelmani [1] core collapse simulation packageis an example of a 3D core-collapse supernova code. It is basedon the open-source
Cactus framework, uses 3D AMR Carte-sian and multi-block grids, and employs many componentsprovided by the open-source
Einstein Toolkit (http://einsteintoolkit.org).
Zelmani has fully general-relativisticgravity and implements general-relativistic MHD. Neutrinosare included either via a rather crude energy-averaged leakagescheme that approximates the overall energetics of neutrinoemission and absorption or via a general-relativistic two-moment M1 radiation-transport solver that has recently beendeployed on first simulations [16].In full radiation-hydrodynamics simulations of the core-collapse supernova problem with 8 levels of AMR,
Zelmani exhibits good strong scaling with hybrid-
OpenMP/MPI to , cores on NSF/NCSA Blue Waters. At larger corecounts, load imbalances due to AMR prolongation and syn-chronization operations begin to dominate the execution time. . D. OTT, MANUSCRIPT FOR COMPUTING IN SCIENCE AND ENGINEERING (CISE), DATED: SEPTEMBER 4, 2018 8 ref. 12 x 6 x2 x Fig. 7. Slices from four semi-global 3D simulations of neutrino-driven convection with parameterized neutrino cooling and heating, carried out in a 45 ◦ wedge. The colormap is the specific entropy; blue colors mark low entropy region, red colors correspond to high entropy. Only the resolution is varied.The wedge marked “ref.” is the reference resolution ( ∆ r = 3 . , ∆ θ = ∆ ϕ = 1 . ◦ ) that corresponds to the resolution of present global 3D detailedradiation-hydrodynamics core-collapse supernova simulations. Note how low resolution favors large flow features and how the turbulence breaks down toprogressively smaller features with increasing resolution. This figure uses simulation results of [21] that includes simulations up to 12 times the referenceresolution that were run on , cores of NSF/NCSA Blue Waters. Rendering by David Radice (Caltech). V. M
ULTI -D D
YNAMICS AND T URBULENCE
Even before the first detailed 2D simulations of neutrino-driven core-collapse supernovae became possible in the mid1990s, it was clear that buoyant convection in the protoneutronstar and in the neutrino-heated region just behind the stalledshock breaks spherical symmetry. Neutrino-driven convectionis due to a negative radial gradient in the specific entropy,making the plasma at smaller radii “lighter” than overlyingplasma. This is a simple consequence of neutrino heating beingstrongest at the base of the heating region. Rayleigh-Taylor-like plumes develop from small perturbations and grow tonon-linear convection. This convection is extremely turbulent,since the physical viscosity in the heating region is vanishinglysmall. Neutrino-driven turbulence is anisotropic on large scales(due to buoyancy), mildly compressible (the flow reachesMach numbers of ∼ . ), and only quasi-stationary, because eventually an explosion develops. Nevertheless, it turns outthat Kolmogorov’s description for isotropic, stationary, in-compressible turbulence works surprisingly well for neutrino-driven turbulence (see Figure 6 for a schematic description ofKolmogorov turbulence and [21]).There is something special about neutrino-driven convectionin core-collapse supernovae: unlike convection in globallyhydrostatic stars, neutrino-driven convection occurs on top ofa downflow of outer core material that has accreted throughthe stalled shock and is headed for the protoneutron star. Theconsequence of this is that there is a competition between( i ) the time it takes for a small perturbation to grow tomacroscopic scale to become buoyant and ( ii ) the time ittakes for it to leave the region that is convectively unstable(the heating region) as it is dragged with the backgroundflow toward the protoneutron star. This means that there are . D. OTT, MANUSCRIPT FOR COMPUTING IN SCIENCE AND ENGINEERING (CISE), DATED: SEPTEMBER 4, 2018 9 log(Wavenumber k)log(Ek) ∝ k -5/3inertial rangeinjectionscale dissipationscale Fig. 6. Schematic view of turbulence: kinetic energy is injected into the flow atlarge scales and cascades through the inertial range via non-linear interactionsof turbulent eddies to small scales (high wavenumbers in the spectral domain)where it dissipates into heat. The scaling of the turbulent kinetic energy withwavenumber in the inertial range is ∝ k − / for Kolmogorov turbulence.This scaling is also found in very high-resolution simulations of neutrino-driven convection [21]. three parameters governing the appearance of neutrino-drivenconvection: the strength of neutrino heating, the initial size ofperturbations entering through the shock, and the downflowrate through the heating region. Because of this, neutrino-driven convection is not a given and simulations find that itdoes not develop in some stars.Even in the absence of neutrino-driven convection, thereis another instability that breaks spherical symmetry in thesupernova core: the standing accretion shock instability (SASI,[3]). SASI was first discovered in simulations that did notinclude neutrino heating. It works via a feedback cycle:small perturbations enter through the shock, flow down tothe protoneutron star and get reflected as sound waves thatin turn perturb the shock. The SASI is a low-mode instabilitythat is most manifest in an up-down sloshing ( (cid:96) = 1 in termsof spherical harmonics) along the symmetry axis in 2D andin a spiral mode ( m = 1 ) in 3D. Once it has reached non-linear amplitudes, the SASI creates secondary shocks (entropyperturbations) and shear flow from which turbulence develops.SASI appears to dominate in situations in which neutrino-driven convection is weak or absent: in conditions whereneutrino heating is weak, the perturbations entering the shockare small, or the downflow rate through the heating region ishigh.Independent of how spherical symmetry is broken in theheating region, all simulations agree that 2D/3D is much morefavorable for explosion than 1D. Some 2D and 3D simulationsyield explosions for stars where 1D simulations fail (see, e.g.,[22]). Why is that?There are two reasons. The first reason has been known forlong and is seemingly trivial: the added degrees of freedom,lateral motion in 2D, and lateral and azimuthal motion in3D, have the consequence that a gas element that entersthrough the shock front spends more time in the heating region before flowing down to settle onto the protoneutronstar. Since it spends more time in the heating region, it canabsorb more neutrino energy, increasing the overall efficiencyof the neutrino mechanism.The second reason has to do with turbulence and hasbecome apparent only in the past few years. Turbulence isoften analyzed employing Reynolds decomposition , a methodthat separates background flow from turbulent fluctuations.Using this method, one can show that turbulent fluctuationslead to an effective dynamical ram pressure (
Reynolds stress )that contributes to the overall momentum balance betweenbehind and in front of the stalled shock. The turbulent pressureis available only in 2D/3D simulations and it has been demon-strated (see, e.g., [23]) that because of this pressure, 2D/3Dcore-collapse supernovae explode with less thermal pressure,and, consequently with less neutrino heating.Now, the Reynolds stress is dominated by turbulent fluctua-tions at the largest physical scales: A simulation that has morekinetic energy in large-scale motions will explode more easilythan a simulation that has less. This realization readily explainsrecent findings by multiple simulation groups: 2D simulationsappear to explode more readily than 3D simulations [22], [23].This is likely a consequence of the different behaviors ofturbulence in 2D and 3D. In 2D, turbulence transports kineticenergy to large scales (which is unphysical), artificially in-creasing the turbulent pressure contribution. In 3D, turbulencecascades energy to small scales (as it should and is knownexperimentally), so a 3D supernova will generally have lessturbulent pressure support than a 2D supernova.Another recent finding by multiple groups is that simula-tions with lower spatial resolution appear to explode morereadily than simulations with higher resolution. There aretwo possible explanations for this and it is likely that theyplay hand-in-hand: ( ) Low resolution creates a numericalbottleneck in the turbulent cascade, artificially trapping tur-bulent kinetic energy at large scales where it can contributemost to the explosion. ( ) Low resolution also increases thesize of numerical perturbations that enter through the shockand from which buoyant eddies form. The larger these seedperturbations are, the stronger is the turbulent convection andthe larger is the Reynolds stress.The qualitative and quantitative behavior of turbulent flow isvery sensitive to numerical resolution. This can be appreciatedby looking at Figure 7, which shows the same 3D simulation ofneutrino-driven convection at 4 different resolutions, spanninga factor of 12 from the reference resolution that is presentlyused in many 3D simulations and which underresolves theturbulent flow. As resolution is increased, turbulent flow breaksdown to progressively smaller features. What also occurs,but cannot be appreciated from a still figure, is that theintermittency of the flow increases as the turbulence is betterresolved. This means that flow features are not persistent, butquickly appear and disappear through non-linear interactionsof turbulent eddies. In this way, the turbulent cascade can betemporarily reversed (this is called backscatter in turbulencejargon), creating large-scale intermittent flow features similarto what is seen at low resolution. The role of intermittencyin neutrino-driven turbulence and its effect on the explosion . D. OTT, MANUSCRIPT FOR COMPUTING IN SCIENCE AND ENGINEERING (CISE), DATED: SEPTEMBER 4, 2018 10 Fig. 8. Visualization of the toroidal magnetic field built up by an inversecascade (large-scale dynamo) from small-scale magnetoturbulence in a mag-netorotational core-collapse supernova. Shown is a ×
70 km octant region with periodic boundaries on the x − z and y − z faces. Regionsof strongest positive and negative magnetic field are marked by light blueand yellowish colors. Dark blue and dark red colors mark regions of weakernegative and positive magnetic field. Based on the NSF/NCSA Blue Waterssimulations of M¨osta et al. mechanism remain to be studied.A key challenge for 3D core-collapse supernova simulationsis to provide sufficient resolution so that kinetic energy cas-cades away from the largest scales at the right rate. Resolutionstudies suggests that this may require between twice to tentimes the resolution of current 3D simulations [21]. A ten-foldincrease in resolution in 3D corresponds to a , timesincrease in the computational cost. An alternative may be todevise an efficient sub-grid model that, if included, providesfor the correct rate of energy transfer to small scales. Workin that direction is still in its infancy in the core-collapsesupernova context.VI. M AKING M AGNETARS :R ESOLVING THE M AGNETOROTATIONAL I NSTABILITY
The magnetorotational mechanism relies on the presence ofan ultra-strong ( ∼ − G ) global, primarily toroidal,magnetic field around the protoneutron star. Such a stronglymagnetized protoneutron star is called a protomagnetar .It has been theorized that the magnetorotational instability(MRI, [7]) could generate a strong local magnetic field that could be transformed into a global field by a dynamo process.While appealing, it was not at all clear that this is whathappens. The physics is fundamentally global and 3D andglobal 3D MHD simulations with sufficient resolution tocapture MRI-driven field growth were impossible to performfor core-collapse supernovae.This changed with the advent of Blue Waters-class petascalesupercomputers and is a testament to how increased computepower and capability systems like Blue Waters facilitate scien-tific discovery. In M¨osta et al. × ×
70 km with uniform resolution. Weperformed four simulations to study the MHD dynamics atresolutions of
500 m ( ∼
200 m ,
100 m , and
50 m ( ∼
20 points per MRI wavelength). Since weemployed uniform resolution and no AMR, the simulationsshowed excellent strong scaling. The
50 m simulation was runon , Blue Waters cores. It consumed roughly
Blue Waters node hours ( ∼
48 million CPU hours).Our simulations with
100 m and
50 m resolution resolvethe MRI and show exponential growth of the magnetic field.This growth saturates at small scales within a few millisecondsand is consistent with what one anticipates on the basis ofanalytical estimates. The MRI drives MHD turbulence that ismost prominent in the layer of greatest rotational shear, justoutside of the protoneutron star core at radii of −
30 km .What we did not anticipate is that in the highest-resolutionsimulation (which resolves the turbulence best), an inverse tur-bulent cascade develops that transports magnetic field energytoward large scales. It acts as a large-scale dynamo that buildsup global, primarily toroidal field, just in the way needed topower a magnetorotational explosion. Figure 8 shows the finaltoroidal magnetic field component in our
50 m simulation after
10 ms of evolution time. Regions of strongest positive andnegative magnetic field are marked by yellowish and light bluecolors, respectively, and are just outside the protoneutron starcore. At the time shown, the magnetic field on large scaleshas not yet reached its saturated state. We expect this to occurafter ∼
50 ms , which could not be simulated.The results of M¨osta et al. suggest that the conditionsnecessary for the magnetorotational mechanism are a genericoutcome of the collapse of rapidly rotating cores. The MRI isa weak field instability and will grow to the needed saturationfield strengths from any small seed magnetic field. The nextstep is to find a way to simulate for longer physical timeand with a larger physical domain. This will be necessaryin order to determine the long-term dynamical impact ofthe generated large-scale magnetic field. Such simulationswill require algorithmic changes to improve parallel scaling,facilitate the efficient use of accelerators, and may require evenlarger and faster machines than Blue Waters. . D. OTT, MANUSCRIPT FOR COMPUTING IN SCIENCE AND ENGINEERING (CISE), DATED: SEPTEMBER 4, 2018 11
VII. C
ONCLUDING R EMARKS
Core-collapse supernova theorists have always been amongthe top group of users of supercomputers. The CDCs and IBMsof the 1960s and 1970s, the vector Crays of the 1970s to1990s, the large parallel scalar architectures of the 2000s, andthe current massively parallel SIMD machines all paved thepath of progress for core-collapse supernova simulations.Today’s 3D simulations are rapidly improving in theirincluded macroscopic and microscopic physics. They arebeginning to answer decades-old questions and are allowingus to formulate new questions. There is still much need forimprovement, which will come at no small price in the post-Moore’s-law era of heterogeneous supercomputers.One important issue that the community must addressis the reproducibility of simulations and the verification ofsimulation codes. It still occurs more often than not thatdifferent codes starting from the same initial conditions andimplementing nominally the same physics arrive at quantita-tively and qualitatively different outcomes. In the mid-2000san extensive comparison of 1D supernova codes took placethat provided results that are still being used as benchmarkstoday [13]. Efforts are now underway that will lead to thedefinition of multi-D benchmarks. In addition to code com-parisons, the increasing availability of open-source simulationcodes and routines for generating input physics (e.g., neutrinointeractions) is furthering reproducibility. Importantly, theseopen-source codes now allow new researchers to enter the fieldwithout the need of spending many years on developing basicsimulation technology that already exists.Core collapse is, in essence, an initial value problem.Current simulations, even those in 3D, start from sphericallysymmetric precollapse conditions from 1D stellar evolutioncodes. However, stars rotate and convection in the layerssurrounding the inert iron care is violently aspherical. Theseasphericities have an impact on the explosion mechanism. Inorder for 3D core-collapse supernova simulations to providerobust and reliable results, the initial conditions must bereliable and robust, and will likely require simulating the finalphases of stellar evolution in 3D [20], which is another multi-D, multi-scale, multi-physics problem.Neutrino quantum-kinetics for including neutrino oscil-lations directly into simulations will be an important, butexceedingly algorithmically and computationally challengingaddition to the simulation physics. Formalisms for doing soare under development and first implementations (in spatially1D) simulations may be available in a few years.A single current top-of-the-line 3D neutrino radiation-hydrodynamics simulation can be carried out to ∼ . − after core bounce at a cost of several tens of millionsof CPU hours; and it still underresolves the neutrino-driventurbulence. What is needed now, are many such simulationsfor studying sensitivity to initial conditions such as rotationand progenitor structure and input physics. These simulationsshould be at higher resolution and carried out for longer sothat the longer-term development of the explosion (or collapseto a black hole) and, for example, neutron star birth kicks canbe reliably simulated. Many longer simulations at higher resolution will requiremuch more compute power than is currently available. Thegood news is that the next generation of petascale systemsand, certainly, exascale machines in the next decade willprovide the necessary FLOPS. The bad news: the radicaland disruptive architectural changes necessary on the routeto exascale will require equally disruptive changes in super-nova simulation codes. Already at petascale, the traditionaldata-parallel, linear/sequential execution model of all presentsupernova codes is the key limiting factor of code performanceand scaling. A central issue is the need to communicate manyboundary points between subdomains for commonly employedhigh-order finite difference and finite volume schemes. Withincreasing parallel process count, communication eventuallydominates over computation in current supernova simulations.Since latencies cannot be hidden, efficiently offloadingdata and tasks to accelerators in heterogeneous systems isdifficult for current supernova codes. The upcoming generationof petascale machines such as DOE’s Summit and Sierra,fully embraces heterogeneity. For exascale machines, powerconsumption will be the driver of computing architecture.Current Blue Waters already draws ∼
10 MW of power andthere is not much upwards flexibility for future machines.Unless there are unforeseen breakthroughs in semiconductortechnology that provide increased single-core performance atorders of magnitude lower power footprint, exascale machineswill likely be all-accelerator with hundreds of millions of slow,highly energy efficient cores.Accessing the compute power of upcoming petascale andfuture exascale machines requires a radical departure fromcurrent code design and major code development efforts.Several supernova groups are exploring new algorithms, nu-merical methods, and parallelization paradigms. Discontinu-ous Galerkin (DG) finite elements (e.g., [25]) have emergedas a promising discretization approach that guarantees highnumerical order while minimizing the amount of subdo-main boundary information that needs to be communicatedbetween processes. In addition, switching to a new, moreflexible parallelization will likely be necessary to preparesupernova codes (and other computational astrophysics codessolving similar equations) for exascale machines. A primecontender being considered by supernova groups is task-basedparallelism, which allows for fine-grained dynamical loadbalancing and asynchronous execution and communication.Frameworks that can become task-based backbones of futuresupernova codes already exist, e.g., Charm++ (http://charm.cs.illinois.edu/research/charm), Legion (http://legion.stanford.edu/overview/), and Uintah (http://uintah.utah.edu/).A
CKNOWLEDGMENTS
I acknowledge helpful conversations with and help fromAdam Burrows, Sean Couch, Steve Drasco, Roland Haas,Kenta Kiuchi, Philipp M¨osta, David Radice, Luke Roberts,Erik Schnetter, Ed Seidel, and Masaru Shibata. I thank theYukawa Institute for Theoretical Physics at Kyoto Universityfor hospitality while writing this article. This work is sup-ported by NSF under award nos. CAREER PHY-1151197 . D. OTT, MANUSCRIPT FOR COMPUTING IN SCIENCE AND ENGINEERING (CISE), DATED: SEPTEMBER 4, 2018 12
EFERENCES[1] C. D. Ott, E. Abdikamalov, P. M¨osta, R. Haas, S. Drasco, E. P.O’Connor, C. Reisswig, C. A. Meakin, and E. Schnetter, “General-relativistic Simulations of Three-dimensional Core-collapse Super-novae,”
Astrophys. J. , vol. 768, 115, May 2013.[2] H. A. Bethe and J. R. Wilson, “Revival of a stalled supernova shock byneutrino heating,”
Astrophys. J. , vol. 295, 14, Aug. 1985.[3] H.-T. Janka, “Explosion Mechanisms of Core-Collapse Supernovae,”
Ann. Rev. Nuc. Par. Sci. , vol. 62, 407, Nov. 2012.[4] P. M¨osta, S. Richers, C. D. Ott, R. Haas, A. L. Piro, K. Boydstun,E. Abdikamalov, C. Reisswig, and E. Schnetter, “MagnetorotationalCore-Collapse Supernovae in Three Dimensions,”
Astrophys. J. Lett. ,vol. 785, L29, Apr. 2014.[5] G. S. Bisnovatyi-Kogan, “The Explosion of a Rotating Star As aSupernova Mechanism.”
Astron. Zh. , vol. 47, 813, Aug. 1970.[6] J. M. LeBlanc and J. R. Wilson, “A Numerical Example of the Collapseof a Rotating Magnetized Star,”
Astrophys. J. , vol. 161, 541, Aug. 1970.[7] S. A. Balbus and J. F. Hawley, “A powerful local shear instability inweakly magnetized disks. I—Linear analysis. II—Nonlinear evolution,”
Astrophys. J. , vol. 376, 214, Jul. 1991.[8] A. Wongwathanarat, H. Janka, and E. M¨uller, “Hydrodynamical NeutronStar Kicks in Three Dimensions,”
Astrophys. J. Lett. , vol. 725, L106,Dec. 2010.[9] C. D. Ott, A. Burrows, L. Dessart, and E. Livne, “Two-DimensionalMultiangle, Multigroup Neutrino Radiation-Hydrodynamic Simulationsof Postbounce Supernova Cores,”
Astrophys. J. , vol. 685, 1069, 2008.[10] E. F. Toro,
Riemann Solvers and Numerical Methods for Fluid Dynam-ics . Berlin: Springer, 1999.[11] T. W. Baumgarte and S. L. Shapiro,
Numerical Relativity: SolvingEinstein’s Equations on the Computer . Cambridge, UK: CambridgeUniversity Press, 2010.[12] T. Kuroda, T. Takiwaki, and K. Kotake, “A New Multi-energy NeutrinoRadiation-Hydrodynamics Code in Full General Relativity and Its Ap-plication to the Gravitational Collapse of Massive Stars,”
Astrophys. J.Supp. Ser. , vol. 222, 20, Feb. 2016.[13] M. Liebend¨orfer, M. Rampp, H.-T. Janka, and A. Mezzacappa, “Super-nova Simulations with Boltzmann Neutrino Transport: A Comparisonof Methods,”
Astrophys. J. , vol. 620, 840, Feb. 2005.[14] K. Sumiyoshi, T. Takiwaki, H. Matsufuru, and S. Yamada, “Multi-dimensional Features of Neutrino Transfer in Core-collapse Super-novae,”
Astrophys. J. Supp. Ser. , vol. 216, 5, Jan. 2015.[15] E. O’Connor and S. M. Couch, “Two Dimensional Core-Collapse Su-pernova Explosions Aided by General Relativity with MultidimensionalNeutrino Transport,” submitted to Astrophys. J.; arXiv:1511.07443 , Nov.2015.[16] L. F. Roberts, C. D. Ott, R. Haas, E. P. O’Connor, P. Diener, andE. Schnetter, “General Relativistic Three-Dimensional Multi-GroupNeutrino Radiation-Hydrodynamics Simulations of Core-Collapse Su-pernovae,”
To appear in Astrophys. J.; arXiv:1604.07848 , Apr. 2016.[17] A. Mirizzi, I. Tamborra, H.-T. Janka, N. Saviano, K. Scholberg, R. Bol-lig, L. H¨udepohl, and S. Chakraborty, “Supernova neutrinos: production,oscillations and detection,”
Nuovo Cimento Rivista Serie , vol. 39, pp. 1–112, 2016.[18] A. Vlasenko, G. M. Fuller, and V. Cirigliano, “Neutrino quantumkinetics,”
Phys. Rev. D. , vol. 89, no. 10, 105004, May 2014.[19] A. W. Steiner, M. Hempel, and T. Fischer, “Core-collapse SupernovaEquations of State Based on Neutron Star Observations,”
Astrophys. J. ,vol. 774, 17, Sep. 2013.[20] S. M. Couch, E. Chatzopoulos, W. D. Arnett, and F. X. Timmes,“The Three-dimensional Evolution to Core Collapse of a Massive Star,”
Astrophys. J. Lett. , vol. 808, L21, Jul. 2015.[21] D. Radice, C. D. Ott, E. Abdikamalov, S. M. Couch, R. Haas, andE. Schnetter, “Neutrino-driven Convection in Core-collapse Supernovae:High-resolution Simulations,”
Astrophys. J. , vol. 820, 76, Mar. 2016. [22] E. J. Lentz, S. W. Bruenn, W. R. Hix, A. Mezzacappa, O. E. B. Messer,E. Endeve, J. M. Blondin, J. A. Harris, P. Marronetti, and K. N. Yakunin,“Three-dimensional Core-collapse Supernova Simulated Using a 15 M (cid:12)
Progenitor,”
Astrophys. J. Lett. , vol. 807, L31, Jul. 2015.[23] S. M. Couch and C. D. Ott, “The Role of Turbulence in Neutrino-driven Core-collapse Supernova Explosions,”
Astrophys. J. , vol. 799, 5,Jan. 2015.[24] P. M¨osta, C. D. Ott, D. Radice, L. F. Roberts, R. Haas, and E. Schnetter,“A Large-Scale Dynamo and Magnetoturbulence in Rapidly RotatingCore-Collapse Supernovae,”
Nature , vol. 528, no. 7582, pp. 376–379,Dec. 2015.[25] J. S. Hesthaven and T. Warburton,