AXISYMMETRIC AB INITIO CORE-COLLAPSE SUPERNOVA SIMULATIONS OF 12–25 M_☉ STARS

Colgate & White (1966) originated the idea that core-collapse supernovae (CCSNe) may be neutrino-driven. Because realistically modeling CCSNe is a computationally complex task of tying together hydrodynamics, neutrino transport, general relativity (GR), nuclear burning, and a realistic equation of state (EoS), and since the many feedbacks and the nonlinearity of the process preclude an analytical treatment, progress in the field has been closely tied to advances in computational resources and code development (Mezzacappa 2005; Kotake et al. 2006; Janka 2012). Two decades ago advances of this sort made possible the first generation of CCSNe simulations in two dimensions (2D; Herant et al. 1994; Burrows et al. 1995; Janka & Müller 1996; Mezzacappa et al. 1998; Fryer & Warren 2002). These simulations, by allowing, among other things, neutrino-driven convection below the shock and continued accretion and shock expansion to coexist, gave explosions. In contrast, continually refined spherically symmetric CCSNe simulations failed to explode except for the lightest progenitors (Rampp & Janka 2000; Mezzacappa et al. 2001; Liebendörfer et al. 2001, 2004; Buras et al. 2003). Further investigation has shown the crucial importance of multidimensional effects, including the standing accretion shock instability (SASI; Blondin et al. 2003), in shaping explosions, and in accounting for the many observed CCSNe explosion asymmetries (see, e.g., Janka 2012).

While the first-generation 2D CCSNe simulations were limited by their use of gray neutrino transport or precomputed spectral neutrino transport, current computing resources have permitted second-generation 2D CCSNe simulations with self-consistent, spectral neutrino transport (Buras et al. 2003, 2006a, 2006b; Burrows et al. 2006, 2007; Bruenn et al. 2006, 2009; Marek & Janka 2009; Suwa et al. 2010; Takiwaki et al. 2012; Müller et al. 2012b, 2012a), essential for accurately computing the neutrino–matter coupling at the heart of the neutrino-driven CCSN mechanism. Realistic multidimensional ab initio (neutrino luminosities and spectra computed from detailed spectral neutrino transport rather than prescribed, parameterized, or incomplete spectral transport) CCSNe simulations are essential for elucidating the details of the neutrino-driven supernova mechanism and fully accounting for supernova observables, including explosion energies, neutrino and gravitational wave signatures, nucleosynthesis, compositional mixing, and neutron star recoils. The Garching group has published such simulations in axisymmetry using the Vertex code, incorporating a complete set of neutrino opacities, velocity-dependent transport terms, and neutrino-energy-bin coupling in the neutrino transport, standard EoSs, and GR corrections to Newtonian self-gravity (Buras et al. 2006a, 2006b; Marek & Janka 2009) or the conformal flatness GR approximation (Müller et al. 2012a, 2012b; see Janka et al. 2012; Müller et al. 2013 for an overview). Their successful models either produce weak explosions or were reported before their explosion energies saturated.

The limited number of published ab initio CCSNe simulations with detailed neutrino transport, the quantitative differences we obtain relative to the Garching group, and the specifics of our findings (explosions more broadly consistent with observations) motivate us to present in this Letter an outline of four simulations being performed with our Chimera code. These simulations span progenitors of mass 12–25 M_☉ with details comparable to Marek & Janka (2009). Here we describe some of the general features already manifested in our simulations and how they support the viability of the neutrino-driven CCSN mechanism.

Chimera is a parallel, multi-physics code built specifically for multidimensional simulation of CCSNe. It is a combination of separate codes for hydrodynamics and gravity; neutrino transport and opacities; and nuclear EoS and reaction network, coupled by a layer that oversees data management, parallelism, I/O, and control. The hydrodynamics are evolved via a dimensionally split, Lagrangian-remap (PPMLR) scheme (Colella & Woodward 1984) as implemented in VH1 (Hawley et al. 2012). Self-gravity is computed by multipole expansion (Müller & Steinmetz 1995) replacing the Newtonian monopole with a GR monopole (Marek et al. 2006, Case A). Neutrino transport is computed in the "ray-by-ray-plus" (RbR+) approximation (Buras et al. 2003), where an independent, spherically symmetric transport solve is computed for each angular "ray" (θ-zone). Neutrinos are advected laterally (in the θ-direction) with the fluid and contribute to the lateral pressure gradient where ρ > 10¹² g cm⁻³. The transport solver is an improved and updated version of the multi-group, flux-limited, diffusion transport solver of Bruenn (1985), enhanced for GR (Bruenn et al. 2001), with an additional geometric flux limiter to prevent the over-rapid transition to free streaming of the standard flux limiter. All O(v/c) observer corrections have been included. We solve for all three flavors of (anti-)neutrinos with four coupled species: ν_e, $\bar{\nu }_e$, ν_μτ = {ν_μ, ν_τ}, $\bar{\nu }_{\mu \tau }=\lbrace \bar{\nu }_{\mu },\bar{\nu }_{\tau }\rbrace$, with 20 energy groups each for α = 4–250 MeV, where α is the lapse function and the comoving-frame group-center energy. The neutrino–matter interactions include emission, absorption, and non-isoenergetic scattering on free nucleons (Reddy et al. 1998) with weak magnetism corrections (Horowitz 2002); emission/absorption (electron capture) on nuclei (Langanke et al. 2003); isoenergetic scattering on nuclei, including ion–ion correlations; non-isoenergetic scattering on electrons and positrons; and pair emission from e⁺e⁻-annihilation (Bruenn 1985) and nucleon–nucleon bremsstrahlung (Hannestad & Raffelt 1998). Unlike the Garching group, we do not include the effective mass corrections at high density or the $\nu _{e}\bar{\nu }_e\rightarrow \nu _{\mu \tau }\bar{\nu }_{\mu \tau }$ process. We utilize the K = 220 MeV incompressibility version of the Lattimer & Swesty (1991) EoS for ρ > 10¹¹ g cm⁻³ and a modified version of the Cooperstein (1985) EoS for ρ < 10¹¹ g cm⁻³ where nuclear statistical equilibrium (NSE) applies. To aid the transition to the XNet (Hix & Thielemann 1999) 14-species α-network (α, ¹²C–⁶⁰Zn) used for the non-NSE regions, we have constructed a 17-species NSE solver to be used in place of the Cooperstein EoS for Y_e > 0.46. An extended version of the Cooperstein electron–photon EoS is used throughout.

During evolution the radial zones are gradually and automatically repositioned during the remap step to track changes in the radial structure. To minimize restrictions on the time step from the Courant limit, we "freeze" the lateral hydrodynamics for the inner eight zones during collapse, and, after prompt convection fades, we expand the laterally frozen region to the inner 8–10 km. In the "frozen" region we set v_θ = 0 and skip the lateral hydrodynamic sweep. The full radial hydrodynamics and neutrino transport are always computed to the center of the simulation for all angular rays. A more extensive description of Chimera is under preparation.

Four non-rotating axisymmetric models (designated B12-WH07, B15-WH07, B20-WH07, and B25-WH07, corresponding to progenitor masses of 12, 15, 20, and 25 M_☉), and identical one-dimensional (1D) models, are initialized without artificial perturbation from the pre-supernova progenitors of Woosley & Heger (2007). A grid of 512 non-equally spaced radial zones covers from the stellar center into the oxygen-rich layers. The 2D models employ 256 uniform-sized angular zones from 0 to π.

All 2D models experience similar evolution during collapse and immediately after bounce. Theta velocities of a few cm s⁻¹ develop by bounce, most likely from numerical roundoff. These grow to at most ∼10 km s⁻¹ and ∼100 km s⁻¹ in the pre-shock and post-shock regions, respectively, at a post-bounce time, t_pb, of 10 ms. At t_pb ∼ 2 ms the shock breaks through the ν_e-sphere, initiating a rapid deleptonization front whose inner edge advects inward with the fluid. An unstable lepton gradient results, which is stabilized by the positive entropy gradient laid down by the shock as it initially strengthens. The initial oscillations of the proto-neutron star (PNS) periodically strengthen and weaken the shock, causing entropy-unstable pockets and transient convective episodes behind the shock. By t_pb ∼ 12 ms an extended unstable entropy gradient has developed behind the shock as it weakens, extending down to ∼60 km and driving a brief episode of convection (Figure 1; 12 ms panel). This episode of "prompt" convection has little appreciable effect on the shock, and by t_pb ∼ 40 ms convective activity has ceased.

Zoom In Zoom Out Reset image size

Up to t_pb ∼ 100 ms, the shock trajectories of corresponding 1D and 2D models track each other closely (Figure 2; lower panel). Remarkably, the shock trajectories for different progenitors track each other rather closely as well. This is a consequence of the compensating factors that determine the radius of the shock during its quasi-stationary accretion phase. Referring to Equation (1) of Janka (2012), these factors (graphed for B12-WH07 and B25-WH07 in Figure 3; upper panel) are the PNS radius and the ν-sphere temperature, whose increase increases the shock radius, and the PNS mass and the mass accretion rate at the shock whose increase has the opposite effect. Despite the variations of these factors, their combinations during the quasi-stationary accretion phase (Figure 3; upper panel, solid lines) are quite similar.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At t_pb ∼ 100 ms several evolutionary developments occur in all 2D models, causing the shock trajectories of the 1D and 2D models to diverge. First, the ratio of τ_adv, the advection timescale through the gain region (region of positive net neutrino heating below the shock), to τ_conv, the Brunt–Väisäla convective growth timescale, begins to exceed the critical value τ_adv/τ_conv ≳ 3 (Foglizzo et al. 2006). Convection in the gain region begins (Figure 1; 90 ms panel), driven by neutrino heating and seeded by the prompt convection episode. Neutrino-driven convection pushes the shock out, enlarging the gain region and increasing τ_adv. Second, the ratio of τ_adv to τ_heat (ratio of the thermal energy E_th to the net neutrino heating rate $\dot{Q}_{\nu _{e}}$ in the gain region) begins to exceed unity. (Using the total energy rather than E_th to compute τ_heat increases by ∼30% the time from bounce for τ_adv/τ_heat to exceed unity.) This signals the beginning of rising entropy in the gain region and an eventual shock revival (Buras et al. 2006a; Thompson et al. 2005). Third, the SASI begins to manifest through low-mode distortions of the shock and oscillating patterns of the radial velocity (evident in the animated version of Figure 1), causing, for example, the gain region to swell along the poles and contract in the equatorial regions (Figure 1; 150 ms panel). This distortion mode increases τ_adv/τ_heat in the expanded regions and decreases it in the contracted regions creating outflow plumes and inflow funnels in the former and latter regions, respectively. Thus, while the shocks of the 1D models slowly contract after t_pb ∼ 100 ms, the average shock radii of the 2D models continue to slowly increase. The expansion accelerates after t_pb ∼ 200 ms (Figure 2).

Going from B12-WH07 to B25-WH07, $\dot{Q}_{\nu _{e}}$ increases ∼2–2.5 fold, E_th increases ∼1.5 fold, and τ_adv decreases ∼1.5 fold, consequently τ_adv/τ_heat ≳ 1 at about the same t_pb for all four models. On the other hand, the decrease in τ_adv noticeably delays the onset of neutrino-driven convection in the more massive models—from t_pb ∼ 60 ms (B12-WH07) to t_pb ∼ 100 ms (B25-WH07). Consequently, it appears that neutrino-driven convection precedes the SASI in B12-WH07 but follows it in the B25-WH07, the SASI having time to saturate in the latter model before the onset of convection (cf. Müller et al. 2012a). B15-WH07 and B20-WH07 are intermediate cases and precedence of convection or the SASI is not clear.

A Legendre decomposition of the shock deformation (computed per Blondin & Mezzacappa 2006) indicates that the growth of the SASI is dominated by the l = 2 (quadrupolar) mode in B12-WH07 and B15-WH07, while the l = 1 (dipolar) mode is dominant in B20-WH07 and B25-WH07. The subsequent rise in the amplitude of shock deformation as the shock begins to accelerate outward is always dominated by the l = 1 mode. The shock deformation oscillates, with a period increasing with the progenitor mass, from 18 ms for B12-WH07 to 30 ms for B25-WH07, though period variations are seen in each model.

The neutrino luminosities L_ν and rms energies 〈E_ν〉_rms for all models follow a similar pattern to that of B12-WH07 (Figure 3; lower panel). Following the ν_e-break-out burst, the luminosities of all neutrino species peak between 100 and 200 ms. The ν_e-, $\bar{\nu }_e$-luminosities, which arise both from the core and from the energy released by accreting matter, exhibit a more pronounced peak during the peak of the mass accretion rate than $L_{\nu _{\mu \tau }}$ and $L_{\bar{\nu }_{\mu \tau }}$, which arise more exclusively from the core. After 200 ms there is a rapid falloff in $L_{\nu _{e}}$ and $L_{\bar{\nu }_e}$ as the shock begins to accelerate outward and the mass accretion rate declines. 〈E_ν〉_rms follows the usual hierarchy, with energy increasing from ν_e to $\bar{\nu }_e$ to ν_μτ to $\bar{\nu }_{\mu \tau }$, the latter three becoming separated by only a few MeV after several hundred ms. The split between $\langle E_{\nu _{\mu \tau }}\rangle _{\rm rms}$ and $\langle E_{\bar{\nu }_{\mu \tau }}\rangle _{\rm rms}$ is due to weak magnetism, which increases (decreases) the opacities of ν_μτ ($\bar{\nu }_{\mu \tau }$). Weak magnetism also causes the $\bar{\nu }_e$-luminosity to exceed the ν_e-luminosity at times after bounce for this model, though the ν_e number luminosity is always dominant. 〈E_ν〉_rms for all neutrino species increases slowly as the ν-spheres contract and heat along with the PNS.

Progenitor masses are correlated with the compactness parameter of O'Connor & Ott (2011; see Table 1), and the luminosities of all neutrino flavors increase with this parameter, nearly doubling from B12-WH07 to B25-WH07 between 50 and 250 ms after bounce. 〈E_ν〉_rms, however, shows only a small increase (a few MeV) with compactness. These trends can be attributed to the dependence of the luminosity on $R^{2}_{\nu \hbox{-}{\rm sph}}T^{4}_{\nu \hbox{-}{\rm sph}}$ and the dependence of the 〈E_ν〉_rms only on T_ν-sph, where R_ν-sph and T_ν-sph are radius and temperature of the ν-sphere, both of which increase with the progenitor compactness during the accretion shock phase (Figure 3, upper panel; note R_ν-sph∝R_PNS).

For all the models, L_ν and 〈E_ν〉_rms are initially isotropic. By t_pb = 60 ms and t_pb = 100 ms for the lower and higher mass models, respectively, anisotropies ∼10% in $L_{\nu _{e},\bar{\nu }_e}$ appear on an l ∼ 10 scale, the scale of the entropy-driven convection appearing for these models at that time. The $L_{\nu _{e}}$ ($L_{\bar{\nu }_e}$) maxima are anti-correlated (correlated) with the convective plumes, the plumes having higher entropy but lower electron fraction (higher positron density). By t_pb = 120 ms anisotropies exist at the ∼30% level in all models on a scale of l = 2–3, with $L_{\nu _{e}}$ and $L_{\bar{\nu }_e}$ maxima becoming more correlated with each other and with the convective plumes and accretion funnels that frequently curve beneath them. By t_pb = 200 ms and thereafter there is a strong l = 1 and l = 2 component in the $L_{\nu _{e},\bar{\nu }_e}$ anisotropy, again correlated with the morphology of the expanding hot plumes and accretion funnels (Figure 1; bottom four panels). $L_{\nu _{\mu \tau },\bar{\nu }_{\mu \tau }}$ remain quite isotropic up to 200 ms after bounce, but follow, with smaller amplitude, the anisotropic pattern of $L_{\nu _{e},\bar{\nu }_e}$ thereafter.

At t_pb ∼ 200 ms several parameters indicate imminent shock revival and the onset of explosion for each model. First, around this time (noted as t_gain in Table 1), the mass residing in the gain region begins to rapidly increase, which increases τ_adv and heating efficiency (Janka 2001; Murphy & Burrows 2008; Marek & Janka 2009; Hanke et al. 2012). Second, the τ_adv/τ_heat rises significantly above unity around this time, signaling shock revival (Buras et al. 2006a; Thompson et al. 2005). Third, near this time, the mean shock radius exceeds 500 km (denoted t₅₀₀ in Table 1).

Table 1. Model Summary Table

Property	Models
	B12-WH07		B15-WH07	B20-WH07	B25-WH07
Progenitor mass, MZAMS (M☉)	12		15	20	25
Progenitor compactness, ξ1.75	0.234		0.598	1.10	1.25
Initial simulation domain radius (km)	30000		20000	20000	20000
Simulation enclosed mass (M☉)	2.08		2.75	3.51	4.57
Radial zone count	512		512	512	512
Angular zone count	256		256	256	256
Collapse time (initiation to bounce; ms)	266		335	462	470
Bounce, central density (1014 g cm−3)	3.13		3.45	3.42	3.29
Bounce, central Ye	0.240		0.249	0.248	0.246
Bounce, shock position (M☉)	0.45		0.46	0.47	0.47
Time to growth in gain region mass, tgain (ms)	256		221	195	192
Time to positive diagnostic energy > 0.01 B (ms)	219		207	192	197
Time to mean shock radius of 500 km, t500 (ms)	236		233	208	212
Post-bounce epoch, tpb (ms)	850	500	500	500	500
Post-bounce diagnostic energy, E+ (B)	0.38	0.33	0.66	0.65	0.70
Post-bounce $E^{+}_{\rm ov}$ (B; see the text)	0.29	0.0.19	0.27	0.03	...
Post-bounce $E^{+}_{\rm ov, rec}$ (B; see the text)	0.32	0.23	0.41	0.08	...
Post-bounce mean shock radius (km)	7245	3429	3293	3632	3783
Post-bounce shock deformation, dshock	0.42	0.31	0.52	0.027	0.31
Post-bounce PNS rest mass (M☉)	1.48	1.48	1.66	1.80	1.88

Download table as:  ASCIITypeset image

A useful parameter characterizing the morphology of the expanding shock is the shock deformation parameter,

$$\begin{equation} d_{\rm shock}= \frac{\max (r_{\rm shock}(\theta)\cos \theta)-\min (r_{\rm shock}(\theta)\cos \theta)}{2\times \max (r_{\rm shock}(\theta)\sin \theta)}-1 \end{equation} \tag{ 1 }$$

(Scheck et al. 2006), where r_shock(θ) is the shock radius as a function of the polar angle, θ. Prolate, oblate, and spherical shocks have positive, negative, and vanishing values of d_shock, respectively. Values of d_shock for the models at t_pb = 500 ms (Table 1) indicate that the shock in three models is significantly prolate while for B20-WH07 it is nearly spherical (see also Figure 2, upper panel).

A measure of the explosion energy frequently used is the diagnostic energy, E⁺ (integral of total energy = kinetic + internal + gravitational, where positive; Buras et al. 2006a; Suwa et al. 2010; Müller et al. 2012b). In Figure 4 we show E⁺ for our models, as well as two other measures of the explosion energy: $E^{+}_{\rm ov}=E^{+}+$ overburden, and $E^{+}_{\rm ov, rec}=E^{+}_{\rm ov}+$ recombination energy. The overburden is the total energy of all negative energy zones on the grid above the innermost positive energy zones plus the total energy of the off-grid material; the recombination energy is the energy that would be released if all ⁴He and neutron–proton pairs recombined to ⁵⁶Ni in the positive total energy zones. We expect that $E^{+}_{\rm ov}$ and $E^{+}_{\rm ov, rec}$ bound the eventual observable explosion energy from below and above, respectively.

Zoom In Zoom Out Reset image size

At 500 ms all three explosion energy measures are positive, and increasing, for B12-WH07, B15-WH07, and B20-WH07, but only E⁺ is positive for B25-WH07 due to its large off-grid binding energy of −0.655 B. For B12-WH07 at 850 ms, the post-shock velocities exceed escape velocity over one-quarter of the shock surface; the accretion rate onto the PNS is essentially zero; and the model is beginning to transit from the simultaneous accretion and outflow phase to the wind phase. The B12-WH07 explosion energies at 850 ms are E⁺ = 0.38 B, $E^{+}_{\rm ov}=0.285$ B, and $E^{+}_{\rm ov, rec}=0.32$ B and appear to be leveling off. These energies are broadly consistent with observations for lower mass progenitors (Smartt 2009).

We have presented results of ab initio axisymmetric CCSN simulations with detailed spectral neutrino transport for a suite of four non-rotating models spanning the mass range 12–25 M_☉, through t_pb = 500 ms for B15-WH07, B20-WH07, and B25-WH07, and t_pb = 850 ms for B12-WH07. At 500 ms, all models show clear indications of developing neutrino-driven explosions, aided initially by strong convective and SASI activity, which imparts pronounced prolate shock deformations in B12-WH07, B15-WH07, and B25-WH07, indicating that large explosion asymmetries are probable (in agreement with observation).

The development of explosions in these models, taking 208–236 ms for the mean shock to reach 500 km, is much slower than in older models—for example, the delayed model of Wilson used by Woosley et al. (1994), for which t₅₀₀ ∼ 150 ms, or the 2D gray neutrino transport models of Herant et al. (1994), for which t₅₀₀ ∼ 75 ms. The explosion development times of the Chimera models reported here are comparable to those of the Garching group for their models that use modern progenitors (t₅₀₀ ∼ 250 ms for model G11 from Müller et al. 2012b and t₅₀₀ ∼ 230 ms for model s27 from Müller et al. 2012a), but significantly shorter than for their models that use older progenitors (t₅₀₀ ∼ 600 ms for model G15 of Müller et al. 2012b). This suggests a significant progenitor dependence, which makes detailed comparisons difficult. The most pronounced difference between our models and those of the Garching group is the growth of the explosion energy. For example, once the explosion develops, E⁺ for both their G15 and s27 models (Janka et al. 2012) grow at a rate of ∼0.5 B s⁻¹, a quarter of the rate in our model B15-WH07. The diagnostic energy in the G11 model saturates at E⁺ ∼ 0.03, one-tenth of the value of our model B12-WH07. Future simulations with both groups using identical progenitors will help elucidate the source of these differences.

Final conclusions regarding the viability of the neutrino-driven CCSN mechanism must wait until different groups obtain similar results from detailed ab initio simulations initiated from the same progenitors and carried out in three dimension. Our simulations thus far support the viability of the neutrino-driven supernova mechanism for the lower progenitor masses and are consistent with observations. These simulations are continuing, and we will publish more complete analyses as they complete.

We thank H.-Th. Janka and B. Müller for a very careful reading of this manuscript and helpful comments. This research was supported by the U.S. Department of Energy Offices of Nuclear Physics and Advanced Scientific Computing Research; the NASA Astrophysics Theory and Fundamental Physics Program (grants NNH08AH71I and NNH11AQ72I); and the National Science Foundation PetaApps Program (grants OCI-0749242, OCI-0749204, and OCI-0749248). This research was also supported by the NSF through TeraGrid resources provided by the National Institute for Computational Sciences under grant number TG-MCA08X010; resources of the National Energy Research Scientific Computing Center, supported by the US DoE Office of Science under Contract No. DE-AC02-05CH11231; and an award of computer time from the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program at the Oak Ridge Leadership Computing Facility, supported by the US DoE Office of Science under Contract No. DE-AC05-00OR22725.