Three-Dimensional Supernova Explosion Simulations of 9-, 10-, 11-, 12-, and 13-M$_{\odot}$ Stars

TL;DR

This study uses advanced 3D simulations to explore supernova explosions in stars of 9 to 13 solar masses, revealing a mass-dependent explodability gap and key factors influencing explosion mechanisms.

Contribution

First 3D simulations of 9-13 solar mass stars showing explodability patterns and identifying critical factors for successful supernova explosions.

Findings

Stars 9-12 M☉ explode easily in 3D

Star 13 M☉ does not explode in simulations

Explosions often have dipolar, asymmetric structures

Abstract

Using the new state-of-the-art core-collapse supernova (CCSN) code F{\sc{ornax}}, we have simulated the three-dimensional dynamical evolution of the cores of 9-, 10-, 11-, 12-, and 13-M$_{\odot}$ stars from the onset of collapse. Stars from 8-M$_{\odot}$ to 13-M$_{\odot}$ constitute roughly 50% of all massive stars, so the explosive potential for this mass range is important to the overall theory of CCSNe. We find that the 9-, 10-, 11-, and 12-M$_{\odot}$ models explode in 3D easily, but that the 13-M$_{\odot}$ model does not. From these findings, and the fact that slightly more massive progenitors seem to explode \citep{vartanyan2019}, we suggest that there is a gap in explodability near 12-M$_{\odot}$ to 14-M$_{\odot}$ for non-rotating progenitor stars. Factors conducive to explosion are turbulence behind the stalled shock, energy transfer due to neutrino-matter absorption and neutrino-matter scattering, many-body corrections to the neutrino-nucleon scattering rate, and the presence of a sharp silicon-oxygen interface in the progenitor. Our 3D exploding models frequently have a dipolar structure, with the two asymmetrical exploding lobes separated by a pinched waist where matter temporarily continues to accrete. This process maintains the driving neutrino luminosty, while partially shunting matter out of the way of the expanding lobes, thereby modestly facilitating explosion. The morphology of all 3D explosions is characterized by multiple bubble structures with a range of low-order harmonic modes. Though much remains to be done in CCSN theory, these and other results in the literature suggest that, at least for these lower-mass progenitors, supernova theory is converging on a credible solution.

Tables

Table 1. Table 1: Some basic model results for the collection of 3D and 2D models calculated for this paper. The model names are followed by the post-bounce time at the end of each simulation, the total explosion energy at the end of each run, the associated “final" baryonic and gravitational masses and proto-neutron-star (PNS) radii (average radius of the 10 11 gm cm -3 surface), and the binding energy of the off-grid stellar envelope. Note that the explosion energies have been corrected for the latter and that after correction all the exploding 3D models have positive explosion energies. The exception is the 13-M ⊙ model, which has yet to explode by the end of both our 2D and 3D simulations.

	t(final)	Exp. Energy	NS mass (bar.)	NS mass (grav.)	PNS Radius	Env. Binding
	(s)	(${10}^{50}$ ergs)	(M⊙)	(M⊙)	(km)	(${10}^{50}$ ergs)
s9.0-2D	1.41	0.71	1.358	1.246	21.3	-0.0206
s9.0-3D	1.042	1.02	1.342	1.233	23.5	$-$
s10.0-2D	1.41	0.62	1.524	1.382	20.8	-0.0953
s10.0-3D	0.767	0.21	1.495	1.358	26.3	$-$
s11.0-2D	1.41	0.78	1.482	1.348	20.8	-0.17
s11.0-3D	0.568	0.75	1.444	1.317	29.7	$-$
s12.0-2D	1.41	1.55	1.568	1.417	20.8	-0.33
s12.0-3D	0.694	0.55	1.507	1.369	28.0	$-$
s13.0-2D	1.311	-0.48	1.854	1.642	21.9	-0.48
s13.0-3D	0.674	-0.48	1.752	1.564	29.1	$-$

Table 2. Table 2: For the runs presented in this paper, the mean shock radius (in units of 1000 kilometers) and mean shock speed (in units of 1000 km s -1 ) at the end of each simulation. Note that the shock is still stalled at the end of the simulation only for the 2D and 3D 13-M ⊙ models.

	t(final)	Mean Shock Radius	Mean Shock Speed
	(s)	(1000 km)	(1000 km s-1)
s9.0-2D	1.41	15.24	14.19
s9.0-3D	1.042	12.42	16.29
s10.0-2D	1.41	7.70	10.62
s10.0-3D	0.767	1.96	6.65
s11.0-2D	1.41	9.18	7.41
s11.0-3D	0.568	2.75	8.00
s12.0-2D	1.41	8.72	8.08
s12.0-3D	0.694	2.66	6.85
s13.0-2D	1.311	0.06	0.067
s13.0-3D	0.674	0.09	0.048

Full text

1 Introduction

Approximately $\sim$50% of the mass function of massive stars above $\sim$8.0 M*⊙* lies below $\sim$13.0 M*⊙. Since only stars more massive than $\sim$8.0 M⊙* can end their lives as core-collapse supernovae (CCSNe) (simultaneously giving birth to either neutron stars or black holes), understanding the mechanism and character of supernova explosions (if they occur) in this modest mass range, assumes an outsized astrophysical importance. Traditionally, those who model the core-collapse and explosion phases of massive stars inherit progenitor models at the cusp of core collapse from experts in massive star evolution. The latter simulate a star’s passage through successive burning phases until an unstable Chandrasekhar white dwarf core emerges at the star’s center, at which point the physical profiles of that core are mapped onto the grid of a supernova code to carry the dynamical, oftimes multi-dimensional, evolution forward. Aside from the stochasticity and chaos associated with the turbulence that attends both progenitor and supernova convective instabilities, the structure of the roughly spherical “initial" model determines the outcome of the supernova simulation, and, it is hoped, the outcome of stellar death. In particular, a progenitor’s radial mass density profile seems to determine much of the subsequent explosive behavior. Figure 1 portrays a representative collection of such profiles. It has been observed, indeed with quantitative variations from modeler to modeler and with some degree of non-monotonic behavior with progenitor ZAMS mass (Woosley & Heger, 2007; Sukhbold et al., 2016; Sukhbold et al., 2018), that massive stars at the lower end of the mass function have steeper mass density profiles with radius than those at the higher end.

It is thought that such steep profiles result in cores that explode easily by the proto-neutron star (PNS) neutrino-driven wind mechanism (Burrows, 1987; Burrows et al., 1995), even in one-dimension (1D, spherical), and this has been shown to be the case (Kitaura et al., 2006; Fischer et al., 2010; Radice et al., 2017) for the pioneering 8.8-M*⊙* model of Nomoto (Nomoto, 1984, 1987). However, such “electron-capture" supernovae (ECSNe) occupy a problematic region of model space (Woosley & Heger, 2015), one in which burning under electron-degenerate conditions with subsequent flashes could 1) eject envelope matter before collapse, 2) radically restructure the core, or 3) compromise the accuracy of 1D stellar evolution simulations. Nevertheless, spherical progenitor models with masses from 8.1 M*⊙* to 9.6 M*⊙,111The 8.1 M⊙* and 9.6 M*⊙* models were for 10*-4* and zero metallicity, respectively. and steep density profiles in the outer Chandrasekhar mantle, have exploded in 1D CCSN simulations (Kitaura et al., 2006; Fischer et al., 2010; Müller et al., 2012; Melson et al., 2015a; Radice et al., 2017). These supernova models universally involve low explosion energies ($\sim$1050 ergs $\equiv$ 0.1 Bethe). When performed in 2D (Burrows et al., 2007; Müller et al., 2012; Radice et al., 2017) or 3D (for the 9.6 M*⊙* model; Melson et al. (2015a)), the explosions are not only low-energy, but quasi-spherical, and such explosions are likely to yield low-mass neutron stars with low-velocity neutron star kicks. A reasonable conclusion is that if the progenitor mass density profiles are as steep as found in these models (Figure 1), the theory and rough explosion numbers arrived at using modern supernova codes that incorporate neutrino transport and heating may be roughly reproducing Nature.

However, the 1D and 2D supernova models of Radice et al. (2017) for the 9, 10, 11 M*⊙* progenitor models of Sukhbold et al. (2016), though spanning a low-mass segment of progenitor parameter space and manifesting a monotonic sequence in mass-density profile from very steep to progressively less steep, do not behave monotonically, nor do they all explode easily. The 9-M*⊙* model explodes easily in 2D by a neutrino-driven wind mechanism (though not in 1D), but the 10-M*⊙* model does not explode in 1D or 2D without significant progenitor velocity perturbations (Couch & Ott, 2013; Müller & Janka, 2015; Müller et al., 2017; Abdikamalov et al., 2016; Takahashi et al., 2016). On the other hand, as with the 9-M*⊙* model, the 11-M*⊙* model does not explode in 1D, but does explode in 2D, and without the aid of perturbations. Moreover, the 12-M*⊙* model of Woosley & Heger (2007), simulated in 2D by Burrows et al. (2018), Vartanyan et al. (2018) and O´Connor & Couch (2018b), does not explode at all unless aided by such things as rotation or significant velocity perturbations (Vartanyan et al., 2018), but the 16-M*⊙* progenitor model from Woosley & Heger (2007) with a significantly shallower mass density profile explodes easily in 2D and 3D (Vartanyan et al., 2019). Collectively, this behavior with variations in mass density profile along the progenitor mass continuum, for which in this mass range the “compactness" parameter (O’Connor & Ott, 2011, 2013; Ott et al., 2018) is monotonic demonstrates once again that the compactness parameter is not predictive of “explodability" (Burrows et al., 2018).

Supernova modeling experience in the recent past now suggests that as the progenitor mass density profile shallows explodability by the wind mechanism is at first easy, and then difficult or impossible. However, with still shallower profiles, explosions by a convection-aided, neutrino-driven mechanism become easier again, though this is not a classic wind. For these progenitors, the post-bounce, pre-explosion mass accretion rates are large enough to result in larger mass densities between the stalled shock wave and PNS core. The correspondingly larger neutrino absorption optical depths in this mantle region lead to more efficient neutrino heating that, with the aid of multi-dimensional effects (Herant et al., 1994; Burrows et al., 1995; Murphy & Burrows, 2008), ignites an explosion. These massive stars have the potential to yield larger explosion energies, with larger kick speeds and neutron star masses. Such explosions are also aided by the accretion of steep silicon/oxygen interfaces (Vartanyan et al., 2018). The associated discontinuity in the mass accretion rate through the shock results in a decrease in the inhibiting tamp, but there is a delay in the corresponding decrease in the neutrino luminosity due to the finite advection time between the shock and the inner core. By accompanying a temporarily unaltered neutrino luminosity with an immediate decrease in the ram pressure, the shock wave can be kicked into explosion by suddenly achieving the critical condition (Burrows & Goshy, 1993). For progenitors with initial mass density profiles intermediate between those of these exploding classes, post-bounce accretion smothers the wind, while failing to provide a large enough neutrino optical depth in the mantle. To date, for these progenitor structures, currently residing near $\sim$12 M*⊙* (Woosley & Heger, 2007; Sukhbold et al., 2016; Sukhbold et al., 2018), 1D and 2D supernova models do not explode in the context of the default physics employed, carried out for as long as one second post-bounce, and without initial rotation nor substantial progenitor seed perturbations. However, we reiterate that the progenitor profiles for massive stars have not theoretically converged and, moreover, that the mapping between ZAMS mass and progenitor structure has not been definitively settled (Woosley & Heger, 2007; Sukhbold et al., 2016; Sukhbold et al., 2018; Meakin et al., 2011; Couch et al., 2015; Müller et al., 2016, 2017).

In this paper, using our new supernova code Fornax (see §2 and Skinner et al. (2019), we simulate in three spatial dimensions (3D) the self-consistent behavior at the lower-mass end (9-, 10-, 11-, 12-, and 13-M*⊙) of a suite of supernova progenitors from Sukhbold et al. (2016). This work is in part a 3D extension of our earlier 1D and 2D study (Radice et al., 2017).222Fornax has also been employed for a variety of 1D and 2D supernova simulations (Wallace et al., 2016; Skinner et al., 2016; Radice et al., 2017; Vartanyan et al., 2018; Burrows et al., 2018), and recently to follow a 3D CCSN explosion for a 16-M⊙* progenitor $\sim$one-second post-bounce (Vartanyan et al., 2019). As stated earlier, these masses span a large fraction of the mass function.

An earlier generation of pioneering 3D simulations (Hanke et al., 2013; Tamborra et al., 2014; Couch & O’Connor, 2014; Melson et al., 2015a; Melson et al., 2015b; Lentz et al., 2015; Müller, 2015; Takiwaki et al., 2016; Roberts et al., 2016; Ott et al., 2018; O’Connor & Couch, 2018a; Summa et al., 2018; Glas et al., 2018a), employing various necessary simplifications, have set the stage for this study, but 3D simulations are still rare. Moreover, some earlier work did not include inelastic scattering nor velocity-dependent transport (Couch & O’Connor, 2014; Ott et al., 2018; Roberts et al., 2016), did not include inelastic scattering on nucleons (Glas et al., 2018a), employed “ray-by-ray+" dimensional reduction (Hanke et al., 2013; Tamborra et al., 2014; Lentz et al., 2015; Müller, 2015; Melson et al., 2015a; Melson et al., 2015b; Takiwaki et al., 2016; Summa et al., 2018), cut out the inner core, or performed the simulations in the inner core in 1D. Depending upon how large this inner 1D core was, the latter procedure could suppress, partially or in full, PNS convection (Dessart et al., 2006; Radice et al., 2017; Glas et al., 2018b). Such convection has been shown to facilitate the explosion of some models (Radice et al., 2017).

Though still expensive, our 3D Fornax runs avoid all these issues, and are state-of-the-art, but still compromise on aspects of the problem. In particular, we employ the two-moment closure M1 method (Vaytet et al., 2011; Skinner et al., 2019), which is not multi-angle.333No such capability with full physics that can follow the evolution for a physically relevant timescale currently exists. In addition, our energy-group number and spatial resolution should be increased,444We use 12 groups per neutrino species and our $r\times\theta\times\phi$ spatial discretization is $678\times 128\times 256$. perhaps significantly, to ensure the results are converged.555However, the resolution dependence of all extant 3D runs remains a universal issue for the CCSN modeling community and our energy and spatial resolutions are better than or comparable to those of most studies in the 3D literature. Furthermore, we use approximate general relativity (GR) to address the enhanced attractive strength of GR and gravitational redshifts, and not full GR (Roberts et al., 2016; Ott et al., 2018). Finally, we lump the $\nu_{\mu}$, $\bar{\nu}_{\mu}$, $\nu_{\tau}$, and $\bar{\nu}_{\tau}$ neutrinos into one species. It has been suggested that this approach can have its limitations (Bollig et al., 2017). Nevertheless, Fornax belongs to a new generation of CCSN simulation codes now emerging with 3D radiation/hydrodynamic capabilities that, though resource-intensive and perforce incorporating various and sundry approximations, finally promise to address the core-collapse supernova phenomenon with full physical fidelity.

In §2, we briefly summarize our computational setup and the microphysics incorporated into Fornax. Then, in §3 we present the central results and derived quantities of our 3D simulations. Following this, in §4 we discuss the characteristics of the explosion debris field and morphology, and then in §5 we recap our most important findings, summarize the current status of supernova theory, and speculate on productive future directions.

2 Methods and Setup

All the 3D and 2D runs performed in this paper used the multi-group radiation/hydrodynamics code Fornax. To date, results using this new capability have been published in Skinner et al. (2016), Wallace et al. (2016), Radice et al. (2017), Vartanyan et al. (2018), Burrows et al. (2018), Skinner et al. (2019), and Vartanyan et al. (2019). The methodologies and equations solved are described in detail in Skinner et al. (2019). The full suite of microphysics and the approach to approximate general relativity (GR) employed are summarized in Vartanyan et al. (2019), in particular in Appendices A and B, respectively, in that paper. We restrict the calculation of gravity to the monopole term, corrected for approximate GR (Marek et al., 2006) and including gravitational redshifts of the neutrino spectra. The comoving-frame transport is a multi-dimensional variant of the two-moment scheme M1 (Vaytet et al., 2011), with analytic closures for the second and third radiation moments (Skinner et al., 2019). Weak magnetism and recoil effects are handled using the formulae in Horowitz (2002) without truncation. Velocity-dependent transport terms to order $(v/c)$ are included, inelastic energy redistribution for neutrino-electron and neutrino-nucleon scattering are handled using the method of Thompson et al. (2003) and Burrows & Thompson (2004), and detailed neutrino-matter scattering and absorption opacities are calculated using the formulations of Burrows et al. (2006), augmented to include the many-body correction to the neutrino-nucleon scattering rates of Horowitz et al. (2017). The latter progressively decreases these neutral-current scattering opacities with increasing density in a manner that approximately captures the associated many-body physics (Burrows & Sawyer, 1998). The effects of this term on core-collapse physics and phenomenology are explored in Burrows et al. (2018) and Vartanyan et al. (2018).

We simulated the multi-D transport of the electron-type neutrinos ($\nu_{e}$), anti-electron-type neutrinos ($\bar{\nu}_{e}$), and “$\nu_{\mu}$"s, where the latter bundles the $\nu_{\mu}$, $\bar{\nu}_{\mu}$, $\nu_{\tau}$, and $\bar{\nu}_{\tau}$ neutrinos. In 3D, we used twelve energy groups for each species, spanning 1 to 300 MeV for the $\nu_{e}$s and 1 to 100 MeV for the others and in 2D we used twenty energy groups over the same ranges.

The calculations were all performed using the same spherical, dendritic grid (Skinner et al., 2019), which deresolves the angular zoning as the center is approached to roughly maintain the same aspect ratio of the spatial zones in the inner $\sim$20 kilometers and to minimize the Courant penalty of spherical convergence. The zoning was 678$\times$128$\times$256 ($r\times\theta\times\phi$) in 3D and 678$\times$256 ($r\times\theta$) in 2D. The inner zone had a $\Delta{r}$ of 0.5 kilometers, and the radial zoning was roughly constant in the inner $\sim$20 km and logarithmic exterior to $\sim$20 kilometers out to an outer spatial boundary of 20,000 kilometers. The 3D runs were begun 10 milliseconds after bounce and were mapped from a corresponding 1D run and continued forward in 3D.

Though 9, 10, and 11 M*⊙* simulations in Radice et al. (2017) employed the non-rotating progenitors from Sukhbold et al. (2016), as we do here for all our models (see Figure 1), the EOS used in that study was that of Lattimer & Swesty (1991) (LS220). Because the SFHo (Steiner et al., 2013) EOS is still consistent with all known nuclear and astrophysical constraints (Tews et al., 2017), whereas the LS220 EOS no longer is, we have opted to use the SFHo EOS in this study. Therefore, we have redone for this paper the comparison 2D simulations for the 9, 10, and 11 M*⊙* progenitors using the SFHo EOS. Also, whereas Radice et al. (2017) investigated the effects of the many-body corrections of Horowitz et al. (2017) (see also Burrows et al. (2018) and Burrows & Sawyer (1998)) to the neutral-current scattering rates off free nucleons, we here include such corrections by default.

A number of authors have touted the potential role of perturbations due to pre-supernova turbulence (Couch & Ott, 2013; Müller & Janka, 2015; Müller et al., 2017; Abdikamalov et al., 2016; Takahashi et al., 2016), preliminary 3D progenitor models are now emerging (Couch et al., 2015; Müller et al., 2017), and Radice et al. (2017) and Vartanyan et al. (2018) have compared various models with and without perturbations. For our 3D runs, we imposed slight $\ell=10$/$m=1$/$n=4$ velocity perturbations on the initial 3D model between 200 and 1000 kilometers, using the scheme of Müller & Janka (2015), as implemented by Radice et al. (2017). However, the magnitude of these perturbations, meant to seed instabilities that might arise, was only 100 km s*-1*. This is to be compared with the radial speed of tens of thousands of km s*-1* at the beginning of the 3D phase. This prescription imposes a minimal degree of perturbation merely to seed turbulence and seems irrelevant to the viability of the explosions we witness. We impose the same perturbations for the comparison 2D simulations, evaluating them in this case at $\phi=0$. Hence, it is highly unlikely that these perturbations played any substantive role in the subsequent growth of turbulence, and were likely dwarfed in effect by the “grid noise."

This paper is meant in part to be a three-dimensional continuation of the 2D and 1D study in Radice et al. (2017). However, that paper explored 9, 10, 11, 8.1, 8.8, and 9.6 M*⊙* models, while we here focus on the 9, 10, 11, 12, and 13 M*⊙* models of Sukhbold et al. (2016). Since the 8.1 and 9.6 M*⊙* progenitors investigated by Radice et al. (2017) were for $10^{-4}$ and zero metallicity and our goals here are to understand the common core-collapse supernova phenomenon, we drop these models from further consideration. In addition, since there is no model in Sukhbold et al. (2016) that corresponds to the 8.8 M*⊙* progenitor of Nomoto (Nomoto, 1984) or that is lower in mass than 9.0 M*⊙*, the latter mass serves as the lower mass limit of this study.

3 Basic Explosion Results

Tables 1 and 2 provide some useful summary numbers associated with our 3D and 2D core-collapse simulations. For Table 1, the explosion energy, the baryon and gravitational masses of the residual proto-neutron star, the average PNS radius (defined as the average radius interior to an isodensity surface at 1011 g cm*-3*), and the envelope binding energy of the off-grid material exterior to 20,000 kilometers (km) are given at t(final), the time post-bounce at the end of each simulation. All the models, except the 13-M*⊙* model, explode. Table 2 lists the mean shock radii and speeds at the end of each simulation. In 3D, the former range from $\sim$15,200 km to $\sim$2,600 km. Since all the exploding models have diagnostic energies (see below) at the end of the calculation in excess of the off-grid overburden energy (see Table 1), these numbers should continue to increase.

Figure 2 depicts the first 0.55 seconds after bounce of the evolution of the mean shock radius for both 3D (thick) and 2D (thin) models.666Given that 3D is the proper context for simulations in Nature, we emphasize that we have carried and calculated the 2D models in this paper merely for important technical comparison. As an aside, we comment that none of these models explodes in 1D. Except for the 13-M*⊙* model, all the 3D and 2D models explode. As indicated on this figure, all the models evolve similarly in the first $\sim$100 milliseconds. This is a consequence of the characteristic timescale in this study for convection and turbulence to kick in and reflects the very small initial seed perturbations imposed (see §2). It is in part the breaking of symmetry due to the onset of turbulence, with its associated dynamics and turbulent stresses, that enables the divergence of the behaviors of the various stalled shocks. The mean shock radii of the 9-M*⊙* and 11-M*⊙* models never really stall or recede and these models explode within $\sim$100 milliseconds of bounce. As with many of our models, the 9-M*⊙* core explodes immediately upon accretion through the partially stalled shock of its sharp silicon/oxygen interface (Vartanyan et al., 2018). The core of the 12-M*⊙* progenitor requires $\sim$40 milliseconds longer to launch its explosion, while the 10-M*⊙* model takes rather much longer to supernova, but is clearly doing so by $\sim$0.3 and $\sim$0.45 seconds in 3D and 2D, respectively. Generally, the 3D models explode slightly earlier than the 2D models, though for the 12-M*⊙* progenitor the 3D and 2D models are launched at roughly the same time. This behavior is consonant with that found by Vartanyan et al. (2019) for the 16-M*⊙* star from Woosley & Heger (2007), who observed that for this progenitor the 3D explosion preceded the 2D explosion by $\sim$50 milliseconds. We note that using their FMD (fully-multi-dimensional) code AENUS-ALCAR, which did not incorporate redistribution by inelastic scattering nor the many-body correction to neutral-current neutrino-nucleon scattering (Horowitz et al., 2017), Glas et al. (2018a) also found that their 3D model exploded before their 2D model, though it exploded at $\sim$300 ms after bounce, later than we find. This difference could easily reflect differences in numerical approach and microphysics. Why the 3D models explode earlier than the 2D models is unclear, but was anticipated in Dolence et al. (2013).777See also Hanke et al. (2012) and Hanke et al. (2013) for an alternate view.

As stated earlier, the 13-M*⊙* model in this progenitor suite does not explode in either 2D or 3D. As Figure 1 demonstrates, this model not only has a shallower density profile (and, therefore higher post-bounce mass accretion rates), but a more muted density jump at its silicon/oxygen interface relative to that of the other “low-mass" progenitors highlighted in this paper. This interface also resides further out in interior mass. These factors clearly have a bearing on the outcome we find. Due to its steep density profile, the 9-M*⊙* model transitions rather quickly to a neutrino-driven wind explosion (Burrows, 1987) and achieves a low asymptotic explosion energy near 0.1 Bethe. The 10-M*⊙, 11-M⊙, and 12-M⊙* models explode similarly weakly (though quickly), due in part to the low absorbing mass (hence, low neutrino optical depth) in the gain region (Bethe & Wilson, 1985) and the low driving $\nu_{e}$ and $\bar{\nu}_{e}$ neutrino luminosities.888Note again that these models also have more prominent Si/O interfaces. Both the low luminosities (which are powered in this early phase mostly by accretion) and low neutrino optical depths are consequences of the steep core mass density profiles of these 9-M*⊙, 10-M⊙, 11-M⊙, and 12-M⊙* models. However, it seems that for these models the positive effect of post-shock turbulence and neutrino heating outweigh the negative effects of the accretion ram and the critical condition for explosion (Burrows & Goshy, 1993; Murphy & Burrows, 2008) can be met. But as the mass density shallows and if the Si/O interface is more subtle, we witness a decreasing tendency to explode, as manifest by our result for the 13-M*⊙* progenitor. However, as our positive result for the corresponding 16-M*⊙* progenitor (Vartanyan et al., 2019) demonstrates, explodability returns despite a shallower, higher density mantle. For such models the neutrino optical depth and accretion-powered luminosities are more significant and, particularly if the Si/O interface jump is not small, “explodability" seems to return. Therefore, we are seeing a gap in explodability, near 13-M*⊙* for the Sukhbold et al. (2016) progenitor models and our Fornax implementation (§2). We saw a similar gap at 12-M*⊙* in Burrows et al. (2018) and Vartanyan et al. (2018) for the 2D models of those studies with a different progenitor suite (Woosley & Heger, 2007), as did O´Connor & Couch (2018b). Hence, it may well be that

1. the lowest-mass progenitors experience low-energy (one to a few$\times$1050 ergs?) explosions that are partially wind-like; 2) higher-mass progenitors experience higher-energy explosions due to higher driving luminosities and neutrino absorption optical depths in the gain region, despite the higher accretion ram tamp; and 3) there is a gap in explodability.999We still suspect that for the shallowest density profiles and highest outer envelope binding energies the latter may be too much to overcome (Burrows et al., 2018). It may be that some flavor of rotation, or new/better physics and numerics, could close this gap, but currently the presence of a gap in explodability, as indicated by our 3D calculations to date, seems plausible. We emphasize that the “compactness parameter" (O’Connor & Ott, 2011) of the 9-M*⊙* through our 16-M*⊙* (Vartanyan et al., 2019) models we have studied to date is monotonic with progenitor mass, but that the outcomes are not.101010For the 9-, 10-, 11-, 12-, and 13-M*⊙* stars of Sukhbold et al. (2016), the compactness parameters, all small, are 3.83$\times 10^{-5}$, 2.64$\times 10^{-4}$, 7.67$\times 10^{-3}$, 2.23$\times 10^{-2}$, and 5.93$\times 10^{-2}$, respectively.

Figure 3 portrays the “diagnostic" explosion energy (Melson et al., 2015a) evolution for our 9-, 10-, 11-, and 12-M*⊙* 3D and 2D models. Since the 13-M*⊙* model does not explode in either 3D or 2D, it does not appear as one of the plotted lines. By diagnostic energy, we mean the sum of the gravitational, kinetic, thermal, and recombination energies of the ejecta. The ejecta are defined as matter that instantaneously appears unbound by the standard Bernoulli condition, and the diagnostic energy does not include the binding energy of the matter exterior to our 20,000-kilometer outer boundary. The latter is given in Table 1 and we see that in all exploding cases the diagnostic energy has overcome the mantle binding term by the end of all our simulations, and is still climbing for most. In fact, the total asymptotic explosion energy (Table 1) is indeed positive for all our exploding models, though those energies are still modest and seem destined to be no more than $\sim$one to a few$\times$1050 ergs. However, for this low-mass massive star subset, this may be realistic (Morozova et al., 2018). The 9-M*⊙* models have total explosion energies (in 3D and 2D) that have almost asymptoted. We note that our 3D 9-M*⊙* model is one of the only state-of-the-art 3D models in the literature to have approached an asymptotic value.

Figure 4 provides the evolution of both the baryon mass (top) and PNS radius (bottom) for both the 3D (thick) and 2D (thin) models. The final values of the residue’s baryon mass (those achieved by the end of each simulation), along with the corresponding gravitational masses, are given in Table 1. We see that the neutron star left behind in the 3D 9-M*⊙* explosion has a baryonic mass of 1.342 M*⊙* and a gravitational mass of 1.246 M*⊙. These values are on the low end of the observed pulsar mass distribution and reflect the lower accretion rates and early explosion times of this model. As suggested in Table 1, such low masses seem to be generic for the steeper density profiles that obtain in this low-mass supernova progenitor mass range. Since the 13-M⊙* models don’t explode and the associated progenitor density profile is the shallowest of the set (implying higher mass accretion rates), the corresponding 3D and 2D PNS baryon masses are large (1.714 [3D] and 1.854 [2D] M*⊙* at simulation end) and are still growing fast at the termination of the runs.

Müller (2015) observes in his long-term study of the 11.2 M*⊙* model of Woosley et al. (2002) that the explosion energy accumulates more unsteadily in 2D than in 3D. In Figure 3 in the early explosion stages, we too see that the 2D explosion energies accumulate more unsteadily than in 3D, but at later times accumulate smoothly. Moreover, we observe for all our models that the residual PNS masses are similar in 2D and 3D. This is in contradistinction to what was found in Müller (2015) and reflects the similar accretion histories in 2D and 3D we witness; what differences there are seem correlated mostly with the (slight) differences in explosion times. An important difference with the Müller (2015) investigation, other than the different code, EOS, and transport methods employed in each study, is that we include as a default in our microphysics suite the Horowitz et al. (2017) many-body correction to the neutrino-nucleon scattering rates. As noted, this facilitates explosion and may mitigate some of the accretional differences between 2D and 3D highlighted in Müller (2015). However, a full exploration of the phenomena called out in Müller (2015) awaits a more exhaustive study that includes longer-term simulations of a broader set of stellar progenitor masses.

The PNS radii are defined as the average radius at a density of 1011 g cm*-3*. As determined in Radice et al. (2017), due to enhanced outward transport and the alteration of inner entropy and Ye profiles, PNS convection increases this radius beyond what would obtain in 1D. The similarity of the 3D and 2D curves in Figure 4 reflects the similar luminosity and core energy loss rates from the interior in 3D and 2D. This fact is reinforced by Figure 5, which demonstrates that the luminosity evolutions in 2D and 3D are nearly the same. We note in passing that the solid-angle-integrated luminosities in 3D evolve more smoothly in time than in 2D, wherein axial sloshing and hydrodynamic asymmetries are more pronounced.

Figure 6 depicts the neutrino heating efficiency, defined as the ratio of the heating rate due to $\nu_{e}$ and $\bar{\nu}_{e}$ absorption in the gain region and the sum of the $\nu_{e}$ and $\bar{\nu}_{e}$ luminosities. This quantity does not include the heating due to inelastic scattering off electrons and nucleons, though this effect is included in the calculations and amounts to a sub-dominant $\sim$5-15% of the total. We see in Figure 6 that the efficiency ranges from a few to $\sim$8% and is similar in 3D and 2D prior to explosion. For the non-exploding 13-M*⊙* model, the 3D and 2D curves do not much differ during the entire simulation, but for the exploding models, since the explosion times differ in 3D and 2D (generally being later in 2D), the associated curves depart from one another after the first of the models explodes. We note that higher efficiencies do not translate into a greater tendency to explode, as witnessed by the behavior of the 13-M*⊙* model. We also note that the neutrino-driven mechanism is not a “1% phenomenon," as frequently considered, but in this context is closer to a “5$-$6%" phenomenon.

Figure 7 shows the evolution of the ejecta mass (in units of M*⊙) versus time after bounce (in seconds, up to a maximum of one second) for both the 3D (thick) and 2D (thin) models. Only the “qualifying" matter on our computational grid interior to the 20,000-kilometer radius of the initial progenitor models is included in this quantity, but it is likely that all the matter exterior to the 20,000-kilometer outer boundary for the exploding models will in fact be ejected. Some matter for the 13-M⊙* models that did not eventually explode still achieved for a short time the Bernoulli condition we have used to identify the ejecta. At the end of the simulations, the ejecta mass (as defined here) for the exploding models ranges from $\sim$10*-2* to $\sim$10*-1* M*⊙*.

Figure 8 renders histograms of the final Ye distributions for the 3D (thick) and 2D (thin) models calculated. This is a true histogram for which the total ejecta mass in a particular Ye bin is the mass given on the y-axis. Importantly, Mej on the ordinate is not a distribution function in Ye (i.e., not $\frac{dM}{dY_{e}}$) and the actual sum of the histogram values given here is the last total ejecta mass plotted in Figure 7 and/or given at the final time listed in Table 1. We find that the total ejecta mass with Ye = 0.5 for the 3D models ranges between $\sim$5$\times$ and $\sim$40$\times$10*-2* M*⊙* and provides a scale for the 56Ni mass ejected. However, we did not perform nucleosynthetic calculations for these runs and are leaving such studies to an upcoming new generation of 3D calculations. As Figure 8 shows, the ejecta for the 3D and 2D runs are not vastly different and are generally proton rich, extending up beyond Ye = 0.56, with most of the ejecta near 0.5. The ejecta Yes are a consequence of the competition between $\nu_{e}$ and $\bar{\nu}_{e}$ absorption and, hence, quite dependent on the fidelity with which the neutrino-matter coupling is handled. This fact should be borne in mind and puts a premium on accurate neutrino transport.

4 Explosion Morphologies

Figures 9 through 12 depict snapshots of volume renderings of the entropy distributions for the 9-M*⊙, 10-M⊙, 11-M⊙, and 12-M⊙* models interior to the shock wave. This post-shock mantle material 1) first participates in the turbulent convection exterior to the PNS core and interior to the stalled shock prior to explosion (left slide) and then 2) constitutes, along with the outer matter encompassed by the expanding shock, the neutrino-driven bubbles and mushroom clouds of the subsequent explosion (right slide). The outer blue shroud is the shock wave. Recall that the 13-M*⊙* model (Figure 13) does not explode and that both its shock wave and PNS core eventually shrink as the core deleptonizes and cools.

Prior to explosion, the characteristic physical scale of the neutrino-driven turbules behind the shock is set by the size of the gain region. This is tens of kilometers to $\sim$100 kilometers and generally represents less than or equal to 50% of the shock radius (to the center). Hence, the angle subtended by the turbules is tens of degrees, which translates into modes of $\ell$ = 5$-$10. To be sure, this is a crude estimate, but serves as a useful mnemonic.111111The growth rates in the earliest phases of instability development will depend upon the initial perturbation spectrum. However, as the explosion commences, the dominant harmonic components of the shock position itself are predominantly the monopole (clearly) and the dipole. In fact, most exploding models have a growing dipole (Dolence et al., 2013), whose growth seems in phase with that of the monopole (whose growth itself constitutes a zeroth-order marker of explosion). For these non-rotating models, the direction chosen by the dipole (a vector) seems arbitrary and is not correlated with the coordinate axes (fortunately), but is correlated with the kick direction (Radice et al. 2019, in preparation). This emphasizes the conclusion that the explosions of such stochastic and chaotic models without physically defined directions would in reality span a uniform spectrum of directions. Figure 14 depicts the monopole-normalized dipole magnitude versus time after bounce for the five 3D models of this paper. The harmonic decomposition employed to generate this figure is that found in Burrows et al. (2012) and Vartanyan et al. (2019). The explosions of our 11-M*⊙* and 12-M*⊙* models are particularly asymmetrical and dipolar, with the dipole component of the shock surface of the 10-M*⊙* model growing fast at the end of that simulation, while by that metric the 3D 9-M*⊙* model is much more spherical. That the latter is more spherical makes sense, since this model explodes almost immediately and a dipole mode doesn’t seem to have time to grow.

Figure 15 renders the corresponding evolution for the normalized higher-order spherical harmonics ($\ell=2,3,4,5$) of the deformation of the shock surface in 3D. For the 3D model results in both Figures 14 and 15, we have lumped together in quadrature the azimuthal $m$s in an effort to avoid the significant clutter that would attend their separate inclusion. However, there is much information in this more refined azimuthal decomposition and we will explore in the future ways to represent it.

We observe that for a given progenitor model the magnitude of the higher-order harmonic coefficients generally decreases as $\ell$ increases and that for the non-exploding 13-M*⊙* model all $\ell$s below six remain of rather low relative amplitude. As both Figures 14 and 15 demonstrate, all the non-monopolar components start small and grow with time as turbulent convection builds. The speed with which this occurs reflects in part the character and magnitude of the seed perturbations, which for all our models are quite small. After explosion onset, for all models the higher-order normalized $\ell$ coefficients grow and then saturate, or grow, peak, slightly diminish, and then saturate (or assume a more secular, long-term drift), with the dipole and quadrupole terms usually predominating in relative magnitude. The time at which the normalized amplitudes of the various harmonics peaks correlates with the explosion time, with the coefficients for the 9-M*⊙* model evolving most quickly through its stages. For this most spherical of our 3D models, the amplitudes of the higher-harmonics remain rather small. For the more-delayed 10-M*⊙* model, the dipole term remains small until explosion is fully underway, and then strengthens (as it does for those models that explode earlier). Before this phase, its normalized shock quadrupole is roughly comparable to its corresponding dipole. Curiously, only the dipole component is of real significance for the 11-M*⊙* model; all its other modes peak early and then assume rather small values. For all models, the temporal fluctuations of all components are significant before explosion, but damp out after explosion as the ejecta evolve into a more coasting, nearly homologous, phase. For the non-exploding 13-M*⊙* model, fluctuations remain rapid on the roughly $\sim$10-millisecond timescales of the turbulence in the gain region. For comparison, we portray in Figure 16 the evolution of the normalized $\ell=1,2,3,4,5$ shock harmonics for the 10-M*⊙* model in 2D. This choice of progenitor is meant merely to be representative. We notice immediately the larger amplitude of the temporal fluctuations in 2D until the simulation is well into explosion, the clear hierarchy in harmonic order with $\ell$, and the larger asymptotic values of the dipole and quadrupole coefficients in 2D than in 3D. The latter are clearly artifacts of the 2D constraint.

We emphasize that all explosions are dominated interior to the shock by bubble structures in mass and entropy, not surprisingly very much like mushroom clouds, with higher harmonic order and angular scales of $\sim$30$-$60 degrees. At late times, but before instabilities in the outer reaches of the massive star can come into play (Utrobin et al., 2018), the bubble structures freeze into a quasi-homologous expansion with frozen angular scales, reflecting the approximate freezing of the associated shock surface structure. However, we note that the influence of computational resolution on the morphology and angular scales of the debris has yet to be determined.

As seen in our 3D 16-M*⊙* explosion model (Vartanyan et al., 2019), two lobes sometimes characterize the debris field, with a “wasp-like" waist separating the two exploding lobes of different sizes (see also Müller (2015), his Figure 6). We see something like this for the 10-M*⊙* and 12-M*⊙* 3D models, but not for the 9-M*⊙* and 11-M*⊙* 3D models and it is curious (perhaps coincidental?) that those former models explode the latest. For a time during the early explosion, this waist, when it appears, is the region where matter is still accreting onto the PNS core in an annular region. This continued accretion helps maintain the driving $\nu_{e}$ and $\bar{\nu}_{e}$ neutrino luminosities, which is emitted more isotropically, despite the directed character of wasp-waist accretion. Such annular accretion seems to funnel matter away from the exploding lobes, thereby diminishing the tamp they experience and, perhaps, facilitating explosion. In some sense, this partial funnelling of matter from antipodal regions into which matter is exploding to an annular region that, at early times, is not exploding is akin to a quadrupolar ($\ell=2$) symmetry breaking in the velocity field and may be a mode found by the outer flow to encourage explosion. Numerous questions emerge, among which are: 1) Is it the case that some models employ such a quadrupolar instability to facilitate explosion? and 2) Given a progenitor, how predictable is the character of the debris field and what is the distribution of explosion morphologies for a given star? We don’t yet know the answers to these questions, but these are worthy future lines of investigation.

The debris fields and morphologies that we have described and that are represented in Figures 9 through 12 are but the initial conditions for the 56Ni-bubble, Rayleigh-Taylor, and Richtmyer-Meshkov instability phases and reverse shocks experienced by the ejecta (Fryxell et al., 1991; Kifonidis et al., 2003; Utrobin et al., 2018) at later times as the blast traverses the outer star. We also note that, aside from the 9-M*⊙* model, our exploding 3D models have not yet asymptoted to their final energies and PNS configurations, though they have achieved positive total explosion energies. Finally, for the non-exploding 13-M*⊙* model, there is some indication at late times of the development of a spiral “SASI" mode (Blondin & Shaw, 2007; Blondin & Mezzacappa, 2007; Rantsiou et al., 2011) and whether this will transition into an explosion at very late times (Takiwaki et al., 2016), though unlikely, has yet to be determined. Other than this late-time manifestation in our 3D 13-M*⊙* run, we don’t witness the classical SASI (Blondin et al., 2003) in any of our simulations. This is in keeping with the conclusion by Burrows et al. (2012) that the SASI is more likely to be in evidence for non-exploding models and that in exploding models neutrino-driven convection generally overwhelms it. However, the final word on the appearance and importance of the SASI has yet to be written.

5 Conclusions

Our new 3D results, and those 1D, 2D, and 3D results of others (Kitaura et al., 2006; Burrows et al., 2007; Fischer et al., 2010; Müller et al., 2012; Melson et al., 2015a; Radice et al., 2017) for this lower-mass range, now suggest that these stars can easily explode by the delayed neutrino mechanism with explosion energies not far from what is observationally expected (Morozova et al., 2018). Moreover, with the results articulated here and elsewhere in the recent literature (e.g., Burrows et al. (2018); Vartanyan et al. (2019)), we can conclude that various features and/or processes disproportionately support or determine explosion in the context of the neutrino heating paradigm. These include 1) turbulence behind the stalled shock wave; 2) the progenitor density structure, in particular the magnitude of the discontinuity at the silicon/oxygen interface; 3) neutrino heating by inelastic neutrino-electron and neutrino-nucleon scattering; 4) many-body corrections to neutrino-nucleon scattering rates (Burrows et al., 2018; Horowitz et al., 2017); and 5) the magnitude and distribution of seed perturbations (Couch & Ott, 2013, 2015; Müller et al., 2017) (though not explored here). We have found in this paper that most low-mass progenitor models explode in 3D, but that there may be a mass gap for the current generation of progenitor models near 12$-$14 M*⊙* where explosion is a bit more problematic. This does not mean that these stars don’t explode, merely that there is a hint from our simulation experience using state-of-the-art tools and modern stellar progenitors that non-rotating cores in this intermediate mass range are less explodable. In addition, we see in our 3D simulations, as with the previous generation of 2D simulations, that simultaneous explosion and accretion (in different solid-angle sectors at a given time) helps maintain explosion in its early phases in some directions by continuing to power the emergent luminosity with an accretion component in other directions. This can’t happen in 1D, but with such a broken symmetry in multi-D it can. Related to this is a wasp-waist pinched structure seen in some debris distributions and the strong dipolar character to many ejecta fields. Moreover, and again (Burrows et al., 2018), we see that compactness parameter in and of itself is not predictive of explosion.

Though the 3D calculations in this new tranche of full-physics supernova models suggest that we and the community of supernova theorists have reached an important milestone in the theory of core-collapse supernovae with the routine generation of sophisticated 3D explosion simulations, there are a number of caveats that bear listing. First, though not outdone by any extant full-physics 3D runs, the spatial and energy-group resolution of our simulations should be improved $-$ a resolution study is in order. Second, we have implemented and fielded with Fornax an approximate GR variant, wherein the correction to the monopolar strength of gravity has been addressed, along with the gravitational redshift of the neutrinos, but bonafide GR (Kuroda et al., 2016; Ott et al., 2018) has been left to later development. Third, we have employed a many-body correction to the axial-vector coupling term in the neutrino-nucleon scattering rates that is still approximate (Horowitz et al., 2017) and have not incorporated the still-uncalculated corresponding correction to the charged-current absorption rates (Burrows & Sawyer, 1999; Roberts et al., 2012; Roberts & Reddy, 2017). Since we have here and elsewhere (Burrows et al., 2018; Vartanyan et al., 2018) found such effects to be important, this remains an important topic for future investigation. Fourth, for this study we have employed the SFHo EOS and though competitive it is not definitive. The dependence of the outcome of stellar collapse on the EOS remains to be determined in full. Fifth, though we perform truly 3D neutrino transfer solving for two angular moments of the specific intensity with a vector neutrino flux, this is not multi-angle transport. The assumed tensor closure form is an ansatz and the higher-order moments (second and third) are provided analytically (Vaytet et al., 2011). Though for pseudo-spherical problems such an approach can be quite accurate (Richers et al., 2017), ultimately full Boltzmann transport in seven dimensions (three space, three momentum space, and time), currently too expensive, will be necessary. In addition, it may well be that aspects of neutrino oscillations must be factored in, raising yet further the ultimate necessary level of computational complexity. Sixth, the mapping between progenitor structure and ZAMS mass is still in flux. Though palpable and enduring progress has been made these last few decades in determining the character of massive star evolution, since the density, seed perturbation, and rotational profiles of the collapsing cores of massive stars are not yet reliably determined in detail, one should be cautious in declaring what a given mass star might do at the end of its life. In addition, we note that core rotation can affect the outcomes, but that rapid rotation that results in neutron-star spin periods less $\sim$50 milliseconds should be rare (Faucher-Giguère & Kaspi, 2006). Such a rotation period is considered “slow" by CCSN theorists. Nevertheless, there is such a densely-packed range in published massive-star structures for the full range of ZAMS masses (see Figure 1) that we find it difficult to conclude that what is now in the literature does not span a realistic range of progenitor structures. However, and importantly, chaos and stochasticity in the turbulent pre-supernova phases will translate into a spread of outcomes in all the astronomical observables (e.g., explosion energies, kick velocities, final neutron star spins and masses, nucleosynthesis, debris spatial distributions). What these spreads are is currently unknown.

Hence, there still remains much to do in CCSN theory that will keep investigators busy for many years. Nevertheless, with the recent emergence of sophisticated codes such as Fornax to generate multiple 3D simulations every year, we have entered a new era in the study of supernova explosions. No longer is the community limited to one or a few expensive 3D runs per year whose individual import is ambiguous. Now, with adequate computational resources many, many full-physics 3D simulations per year are possible. This enables the broader exploration of parameter space and can lead to a more forgiving theoretical environment. The few inevitable mistakes necessary to make real progress in a science are no longer as dire. The last two decades marked the era of commodity 2D simulations that facilitated real progress in parameter exploration and constituted a palpable leap in overall understanding. After $\sim$50 years of development in physics, technique, and computational capabilities, we have now entered the corresponding era in 3D modeling. As a result, with good reason can we expect in the years to come, further (and final?) leaps in insight into this, one of the last remaining fundamental theoretical challenges of stellar astrophysics.

Acknowledgements

The authors acknowledge fundamental contributions to this effort by Josh Dolence and Aaron Skinner. We also acknowledge help with visualization using VisIt from Viktoriya Morozova, fruitful conversations with Hiroki Nagakura and Sean Couch, Evan O’Connor regarding the equation of state, Gabriel Martínez-Pinedo concerning electron capture on heavy nuclei, Tug Sukhbold and Stan Woosley for providing details concerning the initial models, and Todd Thompson regarding inelastic scattering. We acknowledge support from the U.S. Department of Energy Office of Science and the Office of Advanced Scientific Computing Research via the Scientific Discovery through Advanced Computing (SciDAC4) program and Grant DE-SC0018297 (subaward 00009650). In addition, we gratefully acknowledge support from the U.S. NSF under Grants AST-1714267 and PHY-1144374 (the latter via the Max-Planck/Princeton Center (MPPC) for Plasma Physics). DR acknowledges partial support as a Frank and Peggy Taplin Fellow at the Institute for Advanced Study. This overall research project is part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (awards OCI-0725070 and ACI-1238993) and the state of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications. This general project is also part of the “Three-Dimensional Simulations of Core-Collapse Supernovae" PRAC allocation support by the National Science Foundation (under award #OAC-1809073). Moreover, access under the local award #TG-AST170045 to the resource Stampede2 in the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562, was crucial to the completion of this work. Finally, the authors employed computational resources provided by the TIGRESS high performance computer center at Princeton University, which is jointly supported by the Princeton Institute for Computational Science and Engineering (PICSciE) and the Princeton University Office of Information Technology, and acknowledge our continuing allocation at the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the US Department of Energy (DOE) under contract DE-AC03-76SF00098.