THREE-DIMENSIONAL CORE-COLLAPSE SUPERNOVA SIMULATED USING A 15 M_☉ PROGENITOR

That massive stars explode at the end of their lives is well established observationally (Smartt 2009). Numerical simulation of core-collapse supernovae (CCSNe) has been less consistently successful than nature (e.g., Janka 2012; Hix et al. 2014; Mezzacappa et al. 2015). The ultimate source of the neutrino-driven explosion mechanism is the conversion of the gravitational binding energy of the core, collapsed to a proto-neutron star (proto-NS), and of matter accreted onto the proto-NS into neutrinos that heat material behind the shock—reviving it and expelling the stellar envelope as a supernova. This process is decidedly non-spherical. Neutrino heating above the proto-NS drives neutrino-driven convection (νDC; Bethe 1990; Herant et al. 1992, 1994). Also excited are the low-order modes of the standing accretion shock instability (SASI; Blondin et al. 2003). Asphericities in the shock surface channel the continuing accretion into distinct streams. All of these emergent behaviors are manifestly different with imposed axisymmetry (2D) than without (3D); therefore, we should expect 3D modeling to impact the initiation and subsequent development of CCSNe.

Fully capturing the complex behaviors of the CCSN central engine numerically requires spectral neutrino transport, as the neutrinos are not in equilibrium with the fluid and heating is neutrino energy dependent, coupled to self-gravitating hydrodynamics; i.e., spectral neutrino radiation hydrodynamics (νRHD). Few self-consistent, spatially 3D, spectral-νRHD CCSN simulations have been reported. Takiwaki et al. (2012) computed a low-resolution 3D simulation of an 11.2 M_☉ progenitor through the start of explosion. A 2.7° angular resolution follow-up using the same progenitor (Takiwaki et al. 2014) also initiated explosions that subsequently developed more slowly than corresponding 2D simulations. Hanke et al. (2013) computed a non-exploding 2.0° 3D simulation of a 27 M_☉ progenitor, though the same progenitor and physics did explode in 2D. 3D simulations of 11.2 and 20 M_☉ progenitors with the same code and configuration also failed to explode (Tamborra et al. 2014), though Melson et al. (2015a) found a delayed 3D explosion for a 20 M_☉ model with modified opacities. This pattern of delayed and failed explosions in 3D was suggested by Hanke et al. (2012) and later demonstrated explicitly (Couch 2013; Couch & O'Connor 2014) using simulations with parameterized treatments of neutrino heating and cooling, though other simulations showed only small differences or the opposite pattern (Nordhaus et al. 2010; Dolence et al. 2013; Handy et al. 2014; Fernández 2015). Earlier work by Fryer & Warren (2002, 2004) using non-spectral νRHD found explosions in 2D and 3D to be similar.

In this Letter, we report on the first 440 ms of post-bounce evolution for a CCSN initiated from a 15 M_☉ progenitor in 3D and 2D using our multi-dimensional, νRHD code Chimera ⁹ (S. W. Bruenn et al. 2015, in preparation) with modern neutrino interactions and general relativistic corrections to Newtonian self-gravity. We find that the shock revival occurs earlier by ∼100 ms in the 2D simulation relative to 3D and that the growth of the diagnostic explosion energy is similarly accelerated, potentially resulting in stronger explosions in 2D than 3D. This is the first reported 3D explosion at 15 M_☉, a representative mass often used for comparative studies and the first for Chimera. In Section 2, we summarize our methodology and initial conditions. An overview of the simulations is presented in Section 3, with a focus on the differences between the 2D and 3D simulations in Section 4. We discuss our results in context in Section 5, followed by a summary in Section 6.

Initial conditions are taken from the 15 M_☉ pre-supernova progenitor of Woosley & Heger (2007). The inner region (10,700 km; 2.32 M_☉) is remapped onto 540 radial shells on a logarithmic radial grid ($\delta r/r$) modified to track density gradients. Multi-dimensional simulations were initialized from a 1D simulation at 1.3 ms after bounce by applying a 0.1% random density perturbation over radii 10–30 km, mimicking perturbations seen in simulations evolved through bounce in 2D. The angular grid of the 3D simulation (C15-3D) was initialized with a 180-zone (${\rm{\Delta }}\phi =2^\circ$) ϕ grid and a 180-zone θ grid equally spaced in $\mu \equiv \mathrm{cos}\theta$, i.e., equal solid angle. This θ grid widens the pole-adjacent zones (${\rm{\Delta }}{\ell }={R}_{\mathrm{sph}}{\rm{\Delta }}\phi \mathrm{sin}\theta$) and therefore the time step. We evolve in spherical symmetry inside ${R}_{\mathrm{sph}}=6$ km until 45 ms after bounce (when prompt convection fades), thereafter setting ${R}_{\mathrm{sph}}=8$ km. With this grid, the pole-most zone is ≈8.5° wide, resulting in a minimum length and time step ≈4× larger than for a uniform $2^\circ$ θ grid (e.g., Hanke et al. 2013). At 300 ms after bounce, the θ grid was remapped in the 10 θ zones closest to each pole (≈27°) to uniform spacing (${\rm{\Delta }}\theta =2.7^\circ$) and the ϕ sweep at the pole was replaced by a (ϕ) average—yielding similar time steps. The axisymmetric simulation (C15-2D) uses 270 uniform θ zones (${\rm{\Delta }}\theta =2/3^\circ$).

These are the third series of Chimera simulations (Series-C) and are substantially similar to the Series-B simulations (Bruenn et al. 2013, 2014, hereafter B2013 and B2014). A more extensive description of Chimera can be found in Bruenn et al. (2014). The included microphysics are the same as for the Series-B models including the spherical GR terms in the gravity and transport. We solve the multi-group flux-limited diffusion equations in the ray-by-ray approximation for all three flavors of neutrinos and anti-neutrinos with four coupled species: ${\nu }_{e}$, ${\bar{\nu }}_{e}$, ${\nu }_{\mu \tau }=\{{\nu }_{\mu },{\nu }_{\tau }\}$, ${\bar{\nu }}_{\mu \tau }=\{{\bar{\nu }}_{\mu },{\bar{\nu }}_{\tau }\}$, using 20 logarithmically spaced energy groups $\alpha \epsilon =4$–250 MeV, where α is the lapse function and is the comoving-frame group-center energy. The neutrino–matter interactions used are the full set of B2014. We utilize the Lattimer & Swesty (1991) equation of state (EoS; incompressibility K = 220 MeV) for $\rho \gt {10}^{11}$ ${\text{g cm}}^{-3}$ and an enhanced version of the Cooperstein (1985) EoS for $\rho \lt {10}^{11}$ ${\text{g cm}}^{-3}$, and in outer regions, we use a 14-species $\alpha$-network (Hix & Thielemann 1999).

Relative to the Series-B simulations (B2013; B2014), the neutrino transport solver now corrects for frame differences between shock-adjacent zones when computing the flux and flux gradients (S. W. Bruenn et al. 2015, in preparation), permitting spherically symmetric Chimera simulations to track the late shock retreat of the reference simulation in Lentz et al. (2012). This improvement has a modest effect on the shock stalling radius.

All times are given relative to core bounce. The proto-NS is defined as the volume where $\rho \gt {10}^{11}$ ${\text{g cm}}^{-3}$, and the shocked "cavity" is the volume between the proto-NS and the shock.

After remapping from 1D, the multi-D simulations proceed in a similar fashion: convectively unstable regions left behind by the shock progress through the Fe-core trigger prompt convection inside the proto-NS, similar to the axisymmetric Series-B simulations.

Neutrino heating establishes a heating region extending inward from the shock to the gain surface, where net neutrino heating transitions to net cooling. Starting at ≈80 ms for both multi-D simulations, heating at the base of the gain region creates buoyantly unstable conditions, resulting in convective plumes rising against the continuing inflow. Rising plumes begin to affect the shock surface at ≈95 ms for C15-3D and ≈105 ms for C15-2D, as seen by the separation of the minimum and maximum shock radii (Figure 1, dashed lines). Over the next ∼50 ms, both models become completely convective within the shocked cavity. For C15-3D, this results in a flat mean shock radius, ${R}_{\mathrm{shock}}$, that rises gradually up to ≈280 ms. For C15-2D, ${R}_{\mathrm{shock}}$ oscillates and grows faster, indicating earlier shock revival and explosion. The shock for C15-1D, which lacks multi-dimensional flows, reaches a maximum radius of ≈180 km at ≈80 ms and recedes thereafter, which is typical of 1D CCSN simulations.

Zoom In Zoom Out Reset image size

The shock in C15-2D expands rapidly from ≈230 ms onward (Figure 1), with the diagnostic energy¹⁰ ${E}^{+}$ (Figure 2(a)) simultaneously becoming positive. ${E}^{+}$ surpasses 0.01 B by 250 ms and grows rapidly thereafter. For C15-3D, the first evidence of potential explosion begins with an increased growth of ${R}_{\mathrm{shock}}$ at ≈280 ms, accelerating after ≈350 ms, as the largest buoyant plume expands, leading to a small, but growing, ${E}^{+}$.

Zoom In Zoom Out Reset image size

The explosion is clearly more energetic in C15-2D at all times (Figure 2(a)). We evaluate the growth of ${E}^{+}$ over a common period beginning when ${R}_{\mathrm{shock}}$ exceeds 500 km and ending 45 ms later. For C15-3D, ${R}_{\mathrm{shock}}$ passes 500 km at 393 ms when ${E}^{+}$ is 0.034 B, which grows to 0.067 B at 438 ms when ${R}_{\mathrm{shock}}$ is 735 km. For C15-2D, ${R}_{\mathrm{shock}}$ exceeds 500 km at 278 ms when ${E}^{+}$ is 0.041 B, which grows to 0.147 B at 323 ms when ${R}_{\mathrm{shock}}$ reaches 900 km. Over this 45 ms comparison period, the ${E}^{+}$ growth rate is 0.73 ${\text{B s}}^{-1}$ for C15-3D and 2.35 ${\text{B s}}^{-1}$ for C15-2D, and the growth rate of ${R}_{\mathrm{shock}}$ in C15-2D is nearly double that of C15-3D. These indicate a stronger growth of explosion in 2D. Further assessment of ${E}^{+}$ growth requires significantly longer simulations: full saturation of ${E}^{+}$ took ∼1 s to reach ∼1 B for the equivalent 2D model (B15-WH07) of B2014.

The mass in the shocked cavity for C15-2D diverges from that of C15-3D at ≈150 ms (Figure 2(b), solid) and trends strongly upward by ≈250 ms corresponding to rapid growth in ${R}_{\mathrm{shock}}$ and ${E}^{+}$. Similarly, the mass in the gain region (dashed) grows from ≈220 ms onward. For C15-3D, the turnaround in the mass curves is shallower and later and occurs with less mass in the shocked cavity and gain region—all factors that correlate with weaker (or delayed) explosions.

With successful explosions manifest in both 2D and 3D, the leading question is: "Why do axisymmetric models proceed more rapidly (and more forcefully) to explosion than 3D counterparts?"

From 50 to 150 ms, C15-3D shows larger total heating in the gain region (Figure 3 (a)), arising from greater heating efficiency (Figure 3 (c); η _heat, luminosity divided by net heating rate) and neutrino luminosities (Figure 3(b)), though mass accretion at the gain radius (Figure 3(d); ${\dot{M}}_{\mathrm{gain}}$) is similar, resulting in larger ${R}_{\mathrm{shock}}$ than in C15-2D. At ≈150 ms, the ratio of advection and heating time scales (${\tau }_{\mathrm{adv}}/{\tau }_{\mathrm{heat}}$; Figure 3(e)), defined in B2014, grows past unity in both simulations, signaling the potential for thermal runaway (Buras et al. 2006). In C15-2D, ${\tau }_{\mathrm{adv}}/{\tau }_{\mathrm{heat}}$ grows more rapidly, with large excursions driven by the oscillation of the shock along the pole. During this epoch, ${\dot{M}}_{\mathrm{gain}}$ is larger in C15-2D, with large, positive excursions, correlated with favorable increases in the luminosities, ${\eta }_{\mathrm{heat}}$, and ${\tau }_{\mathrm{adv}}/{\tau }_{\mathrm{heat}}$, continuing through C15-2D shock revival (≈250 ms). This favors earlier development of explosion in C15-2D, even though the luminosities remain higher in C15-3D. For C15-2D, luminosities and heating drop with accretion rate after shock revival at ≈250 ms, while both remain noticeably higher in C15-3D. These measures of accretion, luminosity, and heating are generally consistent with the early development of explosion in C15-2D.

Zoom In Zoom Out Reset image size

Recent work has examined the role turbulence can have in supporting the shock leading up to explosion (e.g., Murphy & Meakin 2011; Abdikamalov et al. 2014; Couch & Ott 2015). In Figure 3(f), we show the kinetic energy associated with turbulence in the gain region. Turbulent kinetic energy

$$\begin{eqnarray}&&{E}_{\mathrm{turb}}^{X}={\displaystyle \int }_{\mathrm{gain}}\displaystyle \frac{1}{2}\rho {| {{\boldsymbol{v}}}^{X}| }^{2}\;{dV}\end{eqnarray} \tag{ 1 }$$

is defined for velocities, ${\boldsymbol{v}}$^X, integrated over the gain region. Lateral turbulent energy, ${E}_{\mathrm{turb}}^{\mathrm{lat}}$, is computed by setting radial velocity ${v}_{r}\equiv 0$ in defining ${{\boldsymbol{v}}}^{\mathrm{lat}}$. Anisotropic turbulent energy ${E}_{\mathrm{turb}}^{\mathrm{an}}$, is computed by removing the radial-shell mean $\langle{v}_{r}\rangle$ from ${{\boldsymbol{v}}}^{\mathrm{an}}$, ${v}_{r}^{\mathrm{an}}\equiv {v}_{r}-\langle{v}_{r}\rangle$. The solid lines show the growth of ${E}_{\mathrm{turb}}^{\mathrm{lat}}$, which begins growing at ≈100 ms, the onset of non-radial motions from convection, and continues afterward. ${E}_{\mathrm{turb}}^{\mathrm{an}}$ (dashed lines) is approximately fourfold larger in both simulations prior to shock revival. Both measures are larger for C15-2D than for C15-3D, consistent with Couch & Ott (2015), who posited that stronger turbulent pressure aids the development of explosions in 2D simulations. It is important to note that while ${E}_{\mathrm{turb}}^{}$ measures the kinetic energy of disordered flow, the relevant driver is convection driven by neutrino heating. This is especially important for ${E}_{\mathrm{turb}}^{\mathrm{an}}$ where accretion streams and rising plumes are both deviations from the mean radial flow at large scales.

For multi-dimensional models (including those discussed above), there is a pre-explosion state with convective plumes rising to and distorting the shock. The flow across the accretion shock is diverted by shock geometry toward the local shock minima and then into accretion streams between the plumes. C15-3D, like previously reported 3D simulations, initially shows a large number of small plumes (Figures 4(a), (b)), whereas C15-2D (and most other 2D simulations) shows rapid development toward only a few large plumes (Figure 4(c)). As in B2014, the primary polar plumes in C15-2D oscillate along the symmetry axis, in a manner consistent with the SASI, while neutrino-heated material continues to flow into plumes from the bottom of the gain region, quickly triggering shock revival. In C15-3D, initially small rising plumes are pushed back by ram pressure at the shock and lack the persistence of the larger polar plumes in 2D models (see animated version of Figure 4(a)).

As C15-3D progresses toward explosion, the angular scale of the plumes slowly grows, and, eventually, the largest of these plumes expands continuously and drives shock revival. The growth of the plume scale can be seen at the surface (animation of Figure 4(a)) and deeper (150 km; Figure 5). The growth of the plume scales clearly precedes shock revival at ≈350 ms, when the shock is clearly expanding, and upflow at 150 km is dominated by a single large plume with a few smaller plumes. Parameterized simulations of Fernández (2015) show mergers of plumes, and those of Handy et al. (2014) also found decreasing plume numbers preceding, and after, shock revival. In 2D (including C15-2D), the axis-focused SASI and the tendency for convective vortices to merge ensures that the "large-plume state" is reached earlier (by 150 ms in C15-2D). In 3D, this state takes time to develop and appears to delay the explosion.

Zoom In Zoom Out Reset image size

We have not yet completed a comprehensive analysis of the driver of plume-scale growth, but some observations from our simulation are revealing. As plumes grow in 3D simulations, the associated accretion streams are displaced farther from the central axis of the plume and are therefore less likely to interrupt the inflow of additional buoyant material, providing a "self-shielding" effect that permits further growth. Additionally, the curvature of the shock surface relative to the radial inflow reduces the effective ram pressure normal to the shock surface at locations farther from the plume's axis and thus favors lateral expansion of the plume. Together, these effects could drive the growth of larger plumes, which eventually are able to expand continuously, with material injected from the bottom of the gain region as in 2D simulations.

The effectiveness of this mechanism may be seen in the growth of the plume in the upper left corner of the animated Figure 4(a), which precedes the onset of explosion in 3D. The growth of this plume (also visible in the 200 and 300 ms panels of Figures 4(a), (b)) eventually covers a significant portion of the shocked volume and shock surface solid angle. The lateral expansion of this plume diverts the accretion flow of material, striking it farther and farther from its axis, which contributes to the formation of a strong accretion region on the opposite side of the proto-NS at ≈350 ms, clearly seen in the increased density accreting opposite the largest plume. The growth of the leading plume, and the strong accretion opposite it, give C15-3D a preferred axis (see also Dolence et al. 2013).

The role of dimensionality in CCSN modeling has been of intense interest recently, starting with neutrino "lightbulb" models that parameterize heating with a fixed luminosity. Nordhaus et al. (2010) found that 3D enhanced the potential for explosions, explaining the effect through increased dwell-time for parcels in the neutrino-heating region. Hanke et al. (2012) were unable to confirm that result and for their base 3° simulations found little difference between 2D and 3D. They found that improved angular resolution enhanced explosion in 2D and inhibited explosions in 3D, which they attributed to the differences in the turbulent cascade and the action of the SASI. Couch (2013) found that 3D diminished the potential for explosions, while following the details of the first study carefully. Subsequent studies with other parameterizations have mixed results. Parameterizations necessarily omit physics that might be critical to the nature of the CCSN mechanism; thus, the importance of 2D versus 3D must be made relative to simulations containing all needed physics. Most νRHD simulations, including these, have shown that 3D diminishes the potential for explosion.

The earlier revival of the shock for C15-2D is consistent with the other νRHD simulations in the literature. Axisymmetric simulations of 11.2, 20 M_☉(Tamborra et al. 2013, 2014; Hanke 2014), and 27 M_☉ progenitors (Hanke et al. 2013) produce explosions where their 3D counterparts do not. Melson et al. (2015a) obtained a delayed explosion in 3D by altering neutrino–nucleon scattering for the 20 M_☉ progenitor. In the 3D 20 M_☉ and 27 M_☉ models, they unambiguously demonstrate the spiral SASI mode identified by Blondin & Mezzacappa (2007), whereas C15-3D resembles the neutrino-dominated simulations of Abdikamalov et al. (2014) and Fernández (2015). In their multi-D simulations, the shock recedes in the manner of 1D simulations. Shock contraction in the Hanke et al. (2013) 27 M_☉ model ends when, after accreting a composition interface that drops the (shock) accretion rate from $\approx 0.75$ to ≈0.25 M_☉ s⁻¹, the shock rapidly expands to near its previous peak. Following expansion, the shock revives in their 2D simulation, but contracts in their 3D simulation. Likewise, in the 20 M_☉ simulation, a drop in accretion from ≈0.8 M_☉ s⁻¹ to ≈0.4 M_☉ s⁻¹ at 250 ms expands ${R}_{\mathrm{shock}}$ without explosion, and in the 11.2 M_☉ simulation, a relatively flat accretion of ≈0.2 M_☉ s⁻¹ results in a gradual decline of the 3D shock radius (Tamborra et al. 2014). In Melson et al. (2015a), the same accretion drop in a 20 M_☉ 3D simulation triggers an explosion using enhanced neutrino heating. As in B2014, sudden decreases in accretion at the shock, ${\dot{M}}_{\mathrm{sh}}$, did not trigger shock revival for this progenitor in 2D. For our models at shock revival, ${\dot{M}}_{\mathrm{sh}}$ is ≈0.7 M_☉ s⁻¹ for C15-2D (≈250 ms) and ≈0.5 M_☉ s⁻¹ for C15-3D (≈350 ms) and smoothly declining (Figure 3(d), dashed lines).

In the multi-dimensional models of Takiwaki et al. (2014), ${R}_{\mathrm{shock}}$ rises smoothly beyond maximum 1D ${R}_{\mathrm{shock}}$, with the 2D simulations growing faster than 3D, though these models include less physics than the work of the Garching group discussed above and than our simulations. In their highest-resolution simulations, all of the 2D runs and the longest-run 3D simulation show steepening of the ${R}_{\mathrm{shock}}$ curve, as seen in our models.

In contrast to the above and to our results, Melson et al. (2015b) find a modest increase in explosion energy for 3D versus 2D for a low-mass 9.6 M_☉ Fe-core progenitor, for which they also obtain an explosion in 1D. They attribute this difference to reduced cooling due to increased dissipation of the accretion flow above the gain surface. Given the progenitor differences (1D explosions, plumes not reaching the shock, etc.), these differences relative to our simulations are not surprising.

For the remaining 3D νRHD simulations, including C15-3D, it appears that there is some sensitivity to neutrino heating and the pre-explosion shock history. The 3D simulations of the Garching group (Hanke et al. 2013; Tamborra et al. 2013, 2014; Melson et al. 2015a) all exhibit declining ${R}_{\mathrm{shock}}$ that are only (temporarily) reversed by sharp declines in accretion. Only by modifying the neutrino opacities to account for strange quarks does the ∼10%–20% increase in neutrino heating turn the shock reversal into an explosion. Ostensibly, the Garching code is most similar to Chimera of all multi-dimensional νRHD CCSN codes, though there are some differences: we do not use the ${\nu }_{e}{\bar{\nu }}_{e}\leftrightarrow {\nu }_{\mu \tau }{\bar{\nu }}_{\mu \tau }$ pair-conversion process that increases cooling (Buras et al. 2003); we use the Cooperstein (1985) nuclear EoS in the dynamically important shocked cavity; they use variable Eddington tensor transport along rays while we use flux-limited diffusion; and there are likely differences in numerical dissipation of flows and turbulence from angular and radial resolution and radial grid motion. The effects of these differences can only resolved by code-to-code comparison.

We have computed two identically initialized CCSN simulations, one in axisymmetry (C15-2D) and the other without imposed symmetry (C15-3D), for 440 ms after core bounce using the Chimera CCSN code. The shock revives fairly quickly (≈250 ms after bounce) for C15-2D, following large ${\ell }=1$ oscillations of the shock from the SASI and buoyant plumes along the poles, similar to previous 2D simulations with Chimera (B2013; B2014). Immediately preceding shock revival, C15-2D exhibits more neutrino heating in the gain region and higher accretion at the gain surface than C15-3D. For C15-3D, the shock revival occurs ∼100 ms later, while accretion and convection continue. Though these simulations do not extend long enough to determine final explosion energies, the energies reached in C15-3D, and their growth, make it likely that a lower final explosion energy will be reached in 3D, or it may take longer to reach a similar energy relative to 2D.

The shock is revived in C15-3D through the lead of a single, large-angle plume that results in expansion with a preferred axis, not wholly dissimilar to the axis-imposed structure of 2D simulations. Based on examination of the buoyant plumes and their angular growth in C15-3D, we speculate that the development of fewer large-scale individual plumes may be necessary in 3D to provide the buoyancy needed to overcome the accretion ram pressure and relaunch the shock.

We see evidence for larger turbulent kinetic energy contributing to the earlier shock revival for C15-2D, but further examination is required. Though these are among the best-resolved spectral-νRHD CCSN simulations reported in the literature, we can not be certain the accretion streams and turbulence are adequately resolved. We plan more extensive coverage of these issues, including an examination of resolution, in subsequent publications.

This research was supported by the U.S. Department of Energy Offices of Nuclear Physics and Advanced Scientific Computing Research; the NASA Astrophysics Theory Program (grant NNH11AQ72I); and the National Science Foundation PetaApps Program (grants OCI-0749242, OCI-0749204, and OCI-0749248). P.M. is supported by the National Science Foundation through its employee IR/D program. The opinions and conclusions expressed herein are those of the authors and do not represent the National Science Foundation. This research was also supported by an award of computer time provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program at the Oak Ridge Leadership Computing Facility (OLCF) and at the Argonne Leadership Computing Facility, which are DOE Office of Science User Facilities supported under contracts DE-AC05-00OR22725 and DE-AC02-06CH11357, respectively. Animation by Mike Matheson at OLCF.