Two-dimensional Core-collapse Supernova Explosions Aided by General Relativity with Multidimensional Neutrino Transport

Core-collapse supernovae mark the end stage of stellar evolution for massive stars. These explosions are initiated in the cores of stars with zero-age main-sequence masses above 8–10 M_⊙, where hydrostatic burning of progressively heavier and heavier elements culminates in the formation of a degenerate and inert iron core with a mass of at least ∼1.2 M_⊙. The stability of this iron core against gravitational collapse is due to electron degeneracy pressure. Once the iron core grows beyond the effective Chandrasekhar mass, which can vary depending on the central entropy and electron fraction (Baron & Cooperstein 1990), the core begins to collapse. The collapse is halted when the matter exceeds nuclear density, at which point the repulsive core of the strong nuclear force can supply enough pressure to stabilize the material once again. The inertia of the collapsing inner core causes it to overshoot the new hydrostatic equilibrium, and the rebounding of the core is responsible for launching a shock wave into the still-collapsing outer mantle of the iron core. For a successful core-collapse supernova, this shock will ultimately unbind most of the mass of the star and help spread the nucleosynthetic products of stellar evolution throughout the host galaxy. However, the shock expansion is hampered by two main sinks of energy. First, the dissociation of the heavy nuclei into neutron and protons as they pass through the shock gives an inefficient conversion of the kinetic energy of the infalling material into thermal energy behind the shock, as a modest fraction must go into the internal energy of the matter. Second, the temperature behind the shock and the abundance of free protons now available for electron capture facilitate energy loss via neutrino emission. This also removes thermal energy from behind the shock. Both of these processes reduce the thermal pressure behind the shock and bring it into equilibrium with the ram pressure of the infalling matter. This causes the shock to stall and become an accretion shock. The goal of core-collapse supernova theory over the past ∼50 yr has been to understand the shock revival mechanism. The most ubiquitous mechanism has been the neutrino heating mechanism (Colgate & White 1966; Bethe & Wilson 1985). Neutrino heating is always found in core-collapse simulations. The intense neutrino radiation field streaming away from the proto–neutron star (PNS) quickly falls out of equilibrium with the matter near the neutrinosphere. Outside this radius there are lingering charged-current interactions of electron-type neutrinos and antineutrinos on matter where a net positive transfer of energy between the neutrinos and the matter occurs. If this energy transfer is sufficient, it can lead to shock expansion and a successful explosion. In spherically symmetric simulations of the collapse of typical iron cores, the amount of neutrino heating is not sufficient to drive an explosion (Liebendörfer et al. 2001; Rampp & Janka 2002). Explosions in 2D and 3D simulations with energy-dependent neutrino transport have been reported over the past ∼5–10 yr (Marek & Janka 2009; Suwa et al. 2010; Müller et al. 2012a, 2012b; Takiwaki et al. 2012; Bruenn et al. 2013, 2016; Lentz et al. 2015; Nakamura et al. 2015; Pan et al. 2016; Summa et al. 2016; Suwa et al. 2016; Burrows et al. 2018); however, in this set of successful explosions there are discrepancies among the results from different codes for otherwise very similar initial conditions and input physics.

Many of the earliest simulations of core-collapse supernovae used or were concerned with a general relativistic (GR) approach (Misner & Sharp 1964; Colgate & White 1966; May & White 1966; Wilson 1971; Bruenn 1985; Burrows 1988), and rightly so since neutron stars are sufficiently compact that GR dramatically effects their equilibrium structure. While the importance of including GR gravity in simulations of core-collapse supernovae has always persisted in the literature and is in use in many current and state-of-the-art multidimensional, core-collapse supernova calculations (Marek & Janka 2009; Kuroda et al. 2012; Müller et al. 2012a; Bruenn et al. 2013, 2016; Ott et al. 2013; Lentz et al. 2015; Skinner et al. 2016; Burrows et al. 2018), many modern simulations have used a purely Newtonian approximation for gravity (Takiwaki et al. 2014; Couch & O'Connor 2014; Handy et al. 2014; Couch & Ott 2015; Dolence et al. 2015; Nakamura et al. 2015; Pan et al. 2016; Suwa et al. 2016). Liebendörfer et al. (2001) extensively compared GR gravity and Newtonian gravity in spherical symmetry with a full Boltzmann neutrino transport solver. While their baseline simulations in both prescriptions of gravity fail to explode in 1D, their conclusion is that, overall, GR is helpful for the development of the core-collapse supernova explosion. This conclusion comes out of serendipitous simulations where incorrect nucleon isoenergetic scattering cross sections were used. In these simulations, when GR gravity was used, the simulation predicted an explosion, but when Newtonian gravity was used, the simulations failed to achieve explosions. Bruenn et al. (2001), Buras et al. (2006b), and Lentz et al. (2012) also compared Newtonian and GR simulations in spherical symmetry with energy-dependent neutrino transport. Like Liebendörfer et al. (2001), they observe higher neutrino luminosities and energies but do not extensively study the differences, in part because both simulations still fail to explode owing to the 1D nature of the simulations. Kuroda et al. (2012) examined the difference between special relativistic hydrodynamics and GR hydrodynamics in both 1D and 3D core-collapse simulations using approximate neutrino transport (a gray M1 scheme). However, while the 3D GR simulations show signs of increased susceptibility to explosion, the simulated time was not long enough to observe an explosion. Müller et al. (2012a) also investigated the influence of Newtonian gravity, GR gravity, and GR effective potential gravity in 2D core-collapse simulations using a variable-Eddington-factor, two-moment, energy-dependent neutrino transport scheme and the ray-by-ray+ approximation (where the neutrino transport is done along radial rays assuming spherical symmetry and with minimal coupling between neighboring rays). For the classic 15 M_⊙ model from Woosley & Weaver (1995), Müller et al. (2012a) find a late and asymmetric explosion when using fully GR gravity, but not with Newtonian gravity, GR effective potential gravity, or reduced set neutrino opacities. This was, and still is, the most direct evidence that GR aids in the explosion mechanism. It is important to validate the hypothesis that GR gravity does indeed lead to more favorable conditions for explosion in modern-day simulations, with updated progenitor models, and loosening the assumption of spherically symmetric.

A number of 2D simulations with similar initial conditions have been performed by several different groups using various different treatments of energy-dependent neutrino transport. The initial conditions, chosen by Bruenn et al. (2013), consist of four models from Woosley & Heger (2007) with zero-age main-sequence (ZAMS) masses of 12, 15, 20, and 25 M_⊙. They employ the equation of state (EOS) of Lattimer & Swesty (1991) with incompressibility parameter K₀ = 220 MeV (hereafter LS220). Simulations of these four models using the same EOS were also presented by Summa et al. (2016) and Skinner et al. (2016). Bruenn et al. (2013), Summa et al. (2016), and Skinner et al. (2016) all use Newtonian hydrodynamics with a GR effective potential (Marek et al. 2006, case A). Dolence et al. (2015) also simulate core collapse in these progenitors, but they use the Shen et al. (2011) EOS and assume purely Newtonian gravity. Suwa et al. (2016) and Pan et al. (2016) also simulate at least some of these models, using the LS220 EOS and Newtonian gravity. While no collaborative comparison has been done for these simulations, they do represent the largest set of recent multidimensional core-collapse supernova models performed by multiple simulation codes with similar initial conditions. In Table 1, we summarize the outcome of these simulations (explosion/no explosion), as well as the time of explosion for successful models. We divide the results into simulations with a GR treatment of gravity and those with a Newtonian treatment of gravity. The results are quite disparate. Bruenn et al. (2013) find early and energetic explosions in all four simulated models, while Summa et al. (2016) find explosions in all models, but some explosions are quite delayed (up to ∼500 ms later than Bruenn et al. 2013). Furthermore, Summa et al. (2016) conclude that the final explosion energies will be significantly lower than Bruenn et al. (2013). Skinner et al. (2016) find no explosions in all but one of their models. For the simulations with Newtonian gravity, Dolence et al. (2015) find no explosions, Suwa et al. (2016) find an explosion in one model, the 12 M_⊙ progenitor, and Pan et al. (2016) simulate only s15 and s20 but find an explosion in both models.

In this paper, we perform core-collapse supernova simulations for this same set of progenitors and employ the LS220 EOS (Lattimer & Swesty 1991; O'Connor & Ott 2010). We use a new multidimensional neutrino transport implementation in the FLASH simulation framework (Fryxell et al. 2000; Dubey et al. 2009). We also systematically explore the differences between core-collapse simulations with Newtonian gravity and with a newly implemented GR effective potential. We include our results in Table 1. Our main conclusion is that a GR treatment of gravity is important and required for accurately simulating the neutrino mechanism in core-collapse supernovae. The assumption of Newtonian gravity underestimates the neutrino luminosity and the mean energies of the emitted spectra; consequently, Newtonian gravity results in lower neutrino heating rates in the gain region. Given the large number of simulation groups recently using Newtonian gravity to study core-collapse supernovae, we find it especially important to systematically demonstrate this point. We also compare our results to other core-collapse supernova simulations with comparable initial conditions. Generally, we find that our explosions occur later than those found in the literature. We conclude that these differences are likely due to the differences in the neutrino opacities and interactions used.

The remainder of this paper is organized as follows. In Section 2, we introduce the methods used for computing the GR effective potential and for our neutrino transport scheme. Details and several tests of these methods are presented in the Appendix. In Section 3 we present our main results, including a detailed comparison to publicly available simulation data from Liebendörfer et al. (2005) and O'Connor (2015), multidimensional results assuming Newtonian gravity, and multidimensional results using GR effective potential gravity. We also briefly discuss the impact of resolution and random perturbations on our models. We discuss our Newtonian versus GREP results in Section 4, as well as compare our results with those available in the literature for the same progenitor models and gravity prescriptions. We add four additional models to our GREP simulations and discuss the results in comparison to parameterized 1D models of core collapse. We conclude in Section 5.

Throughout we use units such that G = c = 1.

2.1. Effective General Relativistic Potential

For our implementation of GR gravity, we do not directly solve Einstein's equations; rather, we take an approximate and much simpler approach. We perform our hydrodynamic simulations using standard Newtonian hydrodynamics. However, for the gravitational potential, we replace the Newtonian monopole contribution to the potential with an effective GR potential empirically derived from modified GR structure equations. The formulation we use was presented in Keil (1997) and Rampp & Janka (2002) and later improved by Marek et al. (2006). In spherical symmetry, the Newtonian gravitational potential can be determined via the differential equation

$$\begin{eqnarray}&&\displaystyle \frac{d{\phi }_{\mathrm{NW}}}{{dr}}=\displaystyle \frac{m(r)}{{r}^{2}},\end{eqnarray} \tag{ 1 }$$

where m(r) is the enclosed baryonic mass determined via

$$\begin{eqnarray}&&\displaystyle \frac{{dm}}{{dr}}=4\pi {r}^{2}\rho ,\end{eqnarray} \tag{ 2 }$$

where ρ is the mass density. For Newtonian gravity, the boundary condition on the enclosed mass is m(r = 0) = 0, while the condition on the potential is that at the outer edge of our spherical gravitational potential domain,

$$\begin{eqnarray}&&{\phi }_{\mathrm{NW},\mathrm{bound}}=-{M}_{\mathrm{grid}}/{r}_{\mathrm{outer}}.\end{eqnarray} \tag{ 3 }$$

Marek et al. (2006) modify these equations, using instead

$$\begin{eqnarray}\displaystyle \frac{d{\phi }_{\mathrm{eff}}}{{dr}} & = & \displaystyle \frac{{m}_{\mathrm{TOV}}(r)+4\pi {r}^{3}(P+{P}_{\nu })}{{r}^{2}{{\rm{\Gamma }}}^{2}}\\ & & \times \left[\displaystyle \frac{\rho +\rho \epsilon +P}{\rho }\right],\end{eqnarray} \tag{ 4 }$$

where P, P_ν, and are the matter pressure, neutrino pressure, and specific internal energy, respectively. We take the same boundary condition as Newtonian gravity, ϕ_eff,bound =−M_grid/r_outer. Γ is taken to be

$$\begin{eqnarray}&&{\rm{\Gamma }}={\left[1+{v}^{2}-\displaystyle \frac{2{m}_{\mathrm{TOV}}(r)}{r}\right]}^{1/2},\end{eqnarray} \tag{ 5 }$$

where v² is the square of the average radial velocity and m_TOV(r) is still the enclosed mass, but rather than Equation (2), Marek et al. (2006) ultimately suggest to take the enclosed mass as the solution to the following differential equation:

$$\begin{eqnarray}&&\displaystyle \frac{{{dm}}_{\mathrm{TOV}}}{{dr}}=4\pi {r}^{2}\left(\rho +\rho \epsilon +{E}_{\nu }+\displaystyle \frac{{v}_{i}{F}_{\nu }^{i}}{{\rm{\Gamma }}}\right)\times {\rm{\Gamma }}.\end{eqnarray} \tag{ 6 }$$

Here m_TOV, similar to the Tolman–Oppenheimer–Volkoff (TOV) mass, accounts for internal matter energy and, in this case, the energy density and momentum density of the neutrinos (E_ν and ${F}_{\nu }^{i}$, respectively). The factor of Γ, which is not a part of the standard TOV mass, was empirically found by Marek et al. (2006) to better reproduce actual GR calculations. This is the "case A" GR effective potential of Marek et al. (2006). In the Appendix, we perform an unstable neutron star migration test, in both 1D and 2D, and obtain results consistent with Marek et al. (2006). We have also tested black hole formation in our FLASH implementation using the 40 M_⊙ progenitor model from Woosley & Heger (2007). We find a black hole formation time that is ∼557 ms, or ∼22 ms later than the full GR simulations from O'Connor (2015). Both of these tests give us confidence that our implementation of the GR effective potential is consistent with that of Marek et al. (2006) and close to true GR gravity.

For practical implementation in FLASH, since we use a cylindrical grid, we must spherically average the matter and neutrino fields to obtain ρ, P, P_ν, , v, and E_ν as a function of radius. We take the grid center as the center of our spherical average. We then numerically solve Equations (6) and (5) iteratively, as Γ depends on m_TOV and vice versa. After solving for these terms, we integrate Equation (4) to obtain the spherically symmetric gravitational potential. For the 2D simulations (both Newtonian gravity and effective GR gravity) presented here we only use the monopole term for determining the total gravitational potential.

2.2. Neutrino Radiation Transport

We implement a truly multidimensional, two-moment, energy-dependent, multispecies, neutrino radiation transport scheme with an analytic closure in FLASH. Our implementation follows closely that of O'Connor (2015), which is based on the formalisms of Cardall et al. (2013) and Shibata et al. (2011). In order to incorporate GR aspects into the transport, we follow the assumptions of Rampp & Janka (2002): we retain the g_tt term (here it is the lapse, α) but neglect the g_rr term. As in O'Connor (2015), we distinguish between the laboratory frame, where the neutrino moments are evolved, and the fluid frame, where the energy of each neutrino bin is defined and the neutrino interactions take place. We explicitly couple neighboring energy bins to account for gravitational redshift and velocity gradient terms that transport neutrinos between energy bins. In this work, we ignore inelastic neutrino scattering processes. With these assumptions, the evolution equations for the zeroth (E; neutrino energy density) and first (Fⁱ; neutrino momentum density) energy-dependent neutrino moments in the laboratory frame are, in spherical symmetry,

$$\begin{eqnarray}&&{\partial }_{t}E+\displaystyle \frac{1}{{r}^{2}}{\partial }_{r}[\alpha {r}^{2}{F}^{r}]-{\partial }_{\nu }[\alpha \nu ({L}^{{rr}}{\partial }_{r}{v}_{r}+{F}^{r}{\partial }_{r}\phi )]\\ &&\quad =\,\alpha (W[\eta -{\kappa }_{a}J]-[{\kappa }_{a}+{\kappa }_{s}]{H}^{t}-{F}^{r}{\partial }_{r}\phi ),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\partial }_{t}{F}^{r}+\displaystyle \frac{1}{{r}^{2}}{\partial }_{r}[\alpha {r}^{2}{P}^{{rr}}]-{\partial }_{\nu }[\alpha \nu ({N}^{{rrr}}{\partial }_{r}{v}_{r}+{P}^{{rr}}{\partial }_{r}\phi )]\\ &&\quad =\,\alpha (-[{\kappa }_{s}+{\kappa }_{a}]{H}^{r}+W[\eta -{\kappa }_{a}J]{v}^{r}\\ &&\qquad +\,\displaystyle \frac{{P}^{\theta \theta }+{P}^{\phi \phi }}{r}-E{\partial }_{r}\phi ),\end{eqnarray} \tag{ 8 }$$

and in cylindrical symmetry,

$$\begin{eqnarray}&&{\partial }_{t}E+\displaystyle \frac{1}{r}{\partial }_{r}[\alpha {{rF}}^{r}]+{\partial }_{z}[\alpha {F}^{z}]\\ &&\qquad -\,{\partial }_{\nu }[\alpha \nu ({L}^{{ij}}{\partial }_{i}{v}_{j}+{F}^{i}{\partial }_{i}\phi )]\\ &&\quad =\,\alpha (W[\eta -{\kappa }_{a}J]-[{\kappa }_{a}+{\kappa }_{s}]{H}^{t}-{F}^{i}{\partial }_{i}\phi ),\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\partial }_{t}{F}^{r}+\displaystyle \frac{1}{r}{\partial }_{r}[\alpha {{rP}}^{{rr}}]+{\partial }_{z}[\alpha {P}^{{zr}}]\\ &&\qquad -\,{\partial }_{\nu }[\alpha \nu ({N}^{{rij}}{\partial }_{i}{v}_{j}+{P}^{{ri}}{\partial }_{i}\phi )]\\ &&\quad =\,-\alpha [{\kappa }_{s}+{\kappa }_{a}]{H}^{r}\,+\,W\alpha [\eta -{\kappa }_{a}J]{v}^{r}\\ &&\qquad +\,\alpha \displaystyle \frac{{P}^{\phi \phi }}{r}-\alpha E{\partial }_{r}\phi ,\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{\partial }_{t}{F}^{z}+\displaystyle \frac{1}{r}{\partial }_{r}[\alpha {{rP}}^{{rz}}]+{\partial }_{z}[\alpha {P}^{{zz}}]\\ &&\qquad -\,{\partial }_{\nu }[\alpha \nu ({N}^{{zij}}{\partial }_{i}{v}_{j}+{P}^{{zi}}{\partial }_{i}\phi )]\\ &&\quad =\,\alpha (-[{\kappa }_{s}+{\kappa }_{a}]{H}^{z}+W[\eta -{\kappa }_{a}J]{v}^{z}-E{\partial }_{z}\phi ).\end{eqnarray} \tag{ 11 }$$

In these equations, P^ij is the second moment of the neutrino distribution function, and L^ij and N^ijk are higher-order tensors used for the explicit energy-space flux calculation. $W=1/{[1-{v}_{i}{v}^{i}]}^{1/2}$ is the Lorentz factor, $\alpha =\exp (\phi )$ is the lapse, and ϕ is the gravitational potential from Equation (4). In our Newtonian simulations these GR terms do not appear in the moment evolution equations (i.e., α = 1; ∂_iϕ = 0). Paramters η, κ_a, and κ_s are the neutrino emissivity, absorption opacity, and scattering opacity, respectively, which depend on the local density, temperature, and electron fraction, as well as the neutrino species and energy. J and H^α are the fluid-frame energy density and momentum density, respectively, and must be written in terms of the laboratory frame values so that the equations can be solved implicitly,

$$\begin{eqnarray}&&J={W}^{2}[E-2{F}^{i}{v}_{i}+{v}^{i}{v}^{j}{P}_{{ij}}],\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{H}^{t}={W}^{3}[-(E-{F}^{i}{v}_{i}){v}^{2}+{F}^{i}{v}_{i}-{v}_{i}{v}_{j}{P}^{{ij}}],\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{H}^{i}={W}^{3}[-(E-{F}^{j}{v}_{j}){v}^{i}+{F}^{i}-{v}_{j}{P}^{{ij}}].\end{eqnarray} \tag{ 14 }$$

For the specific details of our implementation, including the methods used to solve for the closure (P^ij, L^ij, and N^ijk), the calculation of the explicit fluxes (both spatial and energy-space), the calculation of the neutrino interaction coefficients (η, κ_a, and κ_s), and the methods used to couple the radiation to the hydrodynamics, we refer the reader to Appendix B.

2.3. FLASH Simulation Setup

3.1. s15WW95—A Detailed Look

We use the 15 M_⊙ core-collapse supernova progenitor profile from Woosley & Weaver (1995) as a test case, denoted throughout as s15WW95. To make comparisons with previous publicly available work (Liebendörfer et al. 2005), we use the Lattimer & Swesty EOS with an incompressibility of 180 MeV (LS180; Lattimer & Swesty 1991). For this test only, following Liebendörfer et al. (2005), we do not include weak magnetism and recoil corrections for the neutrino opacities. We present the results in Figures 1 and 2. For both figures and in all panels, red lines are collapse calculations using Newtonian gravity, while black lines are calculations using our standard GR effective potential, "case A," from Marek et al. (2006). The green lines (thin dot-dashed) use "case R" from Marek et al. (2006) and are presented for comparison purposes with the data extracted from the publicly available Vertex data presented in Liebendörfer et al. (2005). Finally, the blue (thick dot-dashed) line is from a fully GR calculation using GR1D (data are from O'Connor 2015), while the solid blue line and blue dots are data taken from publicly available Agile-Boltztran data from Liebendörfer et al. (2005). For the red and black lines, the thin dashed lines correspond to data taken from GR1D and the thick dashed lines are from FLASH. We note that all of the GR1D simulations (both GREP and full GR) include inelastic scattering of neutrinos on electrons, while FLASH does not. This aspect is especially important during the collapse phase; hence, we start all FLASH simulations from GR1D profiles at ∼15 ms after bounce. For our s15WW95 FLASH simulations the outer extent of our domain is taken to be 10⁹ cm.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 1 shows the fluid-frame neutrino luminosity and rms energy, respectively, of electron neutrinos (top panels), electron antineutrinos (middle panels), and heavy-lepton neutrinos (bottom panels) extracted from the simulations at 500 km. Figure 2 shows the mass in the gain region (top), the total heating rate in the gain region (middle), and heating efficiency (bottom). The gain region is defined as the region within the shock but outside the neutrinosphere, where there is an increasing internal energy of the matter due to neutrino interactions. The heating efficiency is defined as the ratio of the rate of energy deposition into the gain region (${\dot{Q}}_{\mathrm{heat}}^{\mathrm{gain}}$) to the total electron-type neutrino luminosity at the gain radius. Since the latter is difficult to estimate in multidimensional simulations, we estimate it by adding the heating rate to the lab-frame neutrino luminosity extracted from our simulations,

$$\begin{eqnarray}&&{\eta }_{\mathrm{heat}}=\displaystyle \frac{{\dot{Q}}_{\mathrm{heat}}^{\mathrm{gain}}}{{L}_{{\nu }_{e}}^{\mathrm{lab}}+{L}_{{\bar{\nu }}_{e}}^{\mathrm{lab}}+{\dot{Q}}_{\mathrm{heat}}^{\mathrm{gain}}}.\end{eqnarray} \tag{ 15 }$$

For the data from Liebendörfer et al. (2005), we project their fluid-frame luminosity to the lab frame via ${L}_{\nu }^{\mathrm{lab}}\sim {L}_{\nu }^{\mathrm{fluid}}(1+v/c)/(1-v/c)$.

Generally, the deeper gravitational well of the GR simulations (represented by blue, black, and green lines) allows more gravitational binding energy to be released during the accretion phase and results in a higher electron neutrino luminosity and electron antineutrino luminosity, which are predominately fueled by accretion, when compared to the Newtonian simulations (red lines). The deeper gravitational well also leads to higher matter temperatures at the neutrinospheres and therefore higher rms energies. This gives larger heating rates and heating efficiencies. The "case R" effective gravitational potential (labeled as GR1D GREP-R; green curve) overestimates the correction to the gravitational potential and gives even higher luminosities and energies (Marek et al. 2006). Important for our discussion is that our "case R" results (green lines) closely match the "case R" results from the Vertex code (Liebendörfer et al. 2005) This agreement, along with the effective potential test cases in Appendix A, gives us confidence in our effective potential implementation. The improved effective potential presented by Marek et al. (2006; their "case A") comes much closer to reproducing the neutrino luminosities, and matches very closely the rms energies, of the full GR codes Agile-Boltztran and GR1D.

The differences that are seen between the GR1D and GR1D-GREP collapse results are due to the differences between the Newtonian hydrodynamics + GR effective potential gravity and GR hydrodynamics + GR gravity. These differences are also on the order of the differences between the FLASH simulations and the GR1D-GREP simulations, which differ on the hydrodynamics. With the exception of the heavy-lepton neutrino quantities, which are discussed below, the various calculations differ by ∼2%–5% in the neutrino luminosities and rms energies. The differences in the heating-related quantities are somewhat larger, ∼10%, owing to their nonlinear dependence on the neutrino field.

The largest difference between the GR1D and FLASH results is seen in the heavy-lepton sector, which is a direct result of the neglect of inelastic NES. After bounce, the main effect of NES is to downscatter neutrinos (Thompson et al. 2003). Electron-type neutrinos and antineutrinos have charged-current interactions that keep those neutrinos thermodynamically coupled to the matter until low densities, where NES is not as effective. However, heavy-lepton neutrinos, which lack these charged-current interactions, experience a large scattering-dominated region where NES plays an important role. The FLASH simulations, both GREP and NW, have on average 4–6 MeV higher (∼25%) heavy-lepton ν_x rms energy and ∼2 × 10⁵¹ erg s⁻¹ (∼10%) higher luminosities, consistent with the observations of Thompson et al. (2003). However, within the time simulated, the different ν_x luminosities or the additional heating at higher densities supplied by inelastic NES does not have an influence on the electron-type neutrino luminosities or energies, and therefore the heating quantities. For the purposes of neutrino signal prediction, this assumption is quite severe; however, there is little effect on the early core-collapse supernova dynamics.

3.2. 2D Newtonian Results

We perform a total of six 2D core-collapse simulations using Newtonian gravity. For all simulations, we start from profiles taken from GR1D 15 ms after bounce. We perform four simulations of the models introduced in Table 1, s12, s15, s20, and s25 from Woosley & Heger (2007, hereafter WH07). We use our base resolution and the LS220 EOS. As mentioned in Section 2.1, we assume l_max = 0 in our multipole gravity solver. We have tested this assumption by performing a Newtonian transport simulation with model s12 and with l_max = 24. Finally, we test our effective potential solver by performing a Newtonian simulation with the s12 progenitor and where the gravity is computed via the algorithm in Section 2.1, but replacing the differential equations with the Newtonian equivalents. Theoretically, this is precisely the same as the standard monopole term for the gravity, but numerically calculated in a different way. For both of the latter two test simulations we see no differences other than the stochastic behavior of the convection and turbulence.

In Figure 3, we show a subset of results from our Newtonian simulations. We see a very late time (≳1 s) explosion in two models, s20 and s25, and failure in the other models, s12 and s15. Other than the late-time explosion, the simulations of the four models proceed similarly. In both 1D and 2D, in all four models, the shock initially expands out and reaches ∼160–170 km at ∼100 ms after bounce. Only at this time do the 1D and 2D simulations begin to depart from each other, as convection begins to strongly develop behind the shock and provide additional pressure support for the supernova shock (Murphy et al. 2013; Couch & Ott 2015). In 1D and 2D the shocks begin to recede. Compositional interfaces reach the shock front at ∼250, ∼300, and ∼500 ms for models s20, s25, and s12, respectively. At these times, especially in models s25 and s20, we see a short phase of shock expansion. The shock expansion is more prominent in 2D than in 1D. However, this is not sufficient to reenergize the shock in these simulations. At late times the accretion rate has dropped enough to drive an explosion. However, we caution that at these times the baryonic mass of the PNSs is approaching 2.2 M_⊙. The Newtonian gravity approximation is expected to be far from valid in this regime, as the maximum baryonic mass of the LS220 EOS is 2.41 M_⊙.

Zoom In Zoom Out Reset image size

A useful diagnostic of proximity to explosion is the ratio of the advection time through the gain region to the time it takes to unbind the gain material through neutrino heating. A ratio of 1 can be thought of as a runaway condition where matter is unbound before it can be advected into the cooling region. Quantitatively we define this ratio via

$$\begin{eqnarray}&&\displaystyle \frac{{\tau }_{\mathrm{adv}}}{{\tau }_{\mathrm{heat}}}=\displaystyle \frac{{M}_{\mathrm{gain}}/\dot{M}}{| {E}_{\mathrm{gain}}| /{\dot{Q}}_{\mathrm{heat}}}.\end{eqnarray} \tag{ 16 }$$

We define the gain region to consist of any grid zone with a net positive change in the internal energy of the matter due to neutrino interactions. M_gain is the total mass of those grid zones, $\dot{M}$ is the accretion rate measured at 500 km, ${\dot{Q}}_{\mathrm{heat}}$ is the rate of energy exchange between the neutrino and the matter in the gain region, and ${E}_{\mathrm{gain}}={\sum }_{\mathrm{gain}}{m}_{i}[1/2| {v}_{i}{| }^{2}+{\epsilon }_{i}+{\phi }_{i}]$ is the energy of the matter in the gain region. For the latter, v_i is the matter velocity, ϕ_i is the gravitational potential, and _i is the specific internal energy of the matter. For consistency with other works, we take the zero-point defined in the LS220 EOS for _i. In the top right panel of Figure 3 we show this quantity for our Newtonian models. For most of the time, the WH07 models we simulate are well below τ_adv/τ_heat ∼ 1. Only at late times does this ratio approach 1 for the two exploding models.

3.3. 2D Effective GR Results

In our simulations using the GR effective potential we achieve runaway supernova shocks in our core-collapse simulations for three of the four progenitors models: s15, s20, and s25. These explosions occur at significantly earlier times compared to the Newtonian simulations. We attribute this to the increased heating in the GREP simulations discussed in the context of 1D models presented in Section 3.1. This increased heating is a result of the more compact PNS configurations afforded by GR gravity. This gives higher electron neutrino and antineutrino luminosities and rms energies. In Figure 4, we show the evolution of the mean shock radius in our GREP simulations. The top panel shows the early shock development, while the middle panel shows the shock evolution out to 2000 km. The WH07 simulations start very similarly to the pure Newtonian simulations. Initially, the shock radius extends out to ∼145 km at ∼100 ms, roughly 15 km lower than the Newtonian simulations. As with the Newtonian simulations, the shock begins to recede in both 1D and 2D. By the end of the simulated time, we achieve explosions in three of the four models. The times at which the mean shock radii reach 500 km are ∼1100, ∼820, and ∼670 ms for models s15, s20, and s25, respectively. These values are listed in Table 1. The s15WW95_LS180 model is markedly different from the WH07 models; the shock radius stagnates at a slightly lower radius, ∼140 km, and then extends to a much larger radius (∼170 km) after the compositional interface accretes through the shock at ∼150 ms, before ultimately receding and failing to be reenergized.

Zoom In Zoom Out Reset image size

In the bottom panel of Figure 4, we show the ratio of the advection time through the gain region to the time it takes to unbind the matter in the gain region due to neutrino heating, τ_adv/τ_heat, defined in Equation (16). In models s15, s20, and s25, τ_adv/τ_heat surpasses unity at ∼1 s, ∼700 ms, and ∼560 ms, which corresponds to the time when the mean shock radius is starting its outward propagation. The s12 and s15WW95_LS180 models never surpass the value of unity; these models do not explode. For the exploding models, large, late-time asymmetries in the explosion morphology lead to values of τ_adv/τ_heat less than unity. This is due to the breakdown of the underlying assumptions of Equation (16), which is technically only valid for spherical flows.

The left column of Figure 5 shows the luminosity of each of the three neutrino species. Prior to explosion, or for the entire evolution for models s12 and s15WW95_LS180, the neutrino luminosities closely follow the 1D simulation. This is because the neutrino luminosity being emitted in simulations with stalled or failing supernova shocks is set by the accretion rate, which is the same in both 1D and 2D. The situation changes after the explosion sets in. After this time, not all accreted matter is able to settle into the cooling region since it is being heated by neutrinos and unbound. As a result, the released gravitational binding energy is not being converted into neutrino luminosity. The consequence is that the post-explosion neutrino luminosity in 2D is lower than the neutrino luminosity determined in the 1D simulation at the same time. Our ν_x luminosities are not as dramatically influenced by the fluid motions in the gain region, as they are emitted from deeper in the cooling region away from the convective motions. Consequently, they show less variability than both ν_e and ${\bar{\nu }}_{e}$. However, there is a significant multidimensional effect on the ν_x neutrinos. We see a boost in the ν_x luminosities due to PNS convection. This convection slows the recession of the neutrinosphere and dredges up heat from deeper in the PNS. This results in an increase in the ν_x luminosity and ν_x rms energy (right column of Figure 5) starting at ∼200 ms. It is also responsible for lowering the electron-type neutrino and antineutrino rms energies. We give a more detailed discussion of PNS convection in Section 4.

Zoom In Zoom Out Reset image size

In the top panel of Figure 6, we show the early development of the diagnostic explosion energy (Buras et al. 2006b; Bruenn et al. 2013; Melson et al. 2015),

$$\begin{eqnarray}&&{E}_{\exp }^{\mathrm{diag}}={\int }_{{ \mathcal V }}d{ \mathcal V }\rho {[1/2| v{| }^{2}+(\epsilon -{\epsilon }_{0})+\phi ]}_{\gt 0},\end{eqnarray} \tag{ 17 }$$

where the volume integral is over all zones where the integrand is positive. ϕ is the gravitational potential energy, is the specific internal energy, and ₀ is the specific internal energy of matter with the same density and Y_e, but T = 0. This roughly accounts for the energy that will be released into kinetic energy as the matter cools and the free nucleons recombine into alpha particles and later nuclei. To truly determine this energy, one must track the evolution to much longer times and do a proper treatment of nucleosynthesis, something we are not capable of doing with our current infrastructure. We also note that this diagnostic explosion energy does not account for the binding energy of the material outside the shock (i.e., the "overburden"), which can be significant in the higher-mass models, perhaps upward of 1 B (Bruenn et al. 2016). Our explosion energies are quite small compared to the results of Bruenn et al. (2013), and more along the lines predicted from other successful explosions in 2D and 3D (Suwa et al. 2010; Hanke et al. 2012; Müller et al. 2012a; Melson et al. 2015; Pan et al. 2016; Suwa et al. 2016). We do note that the amount of simulated post-explosion time is short. In our simulations, along with the weak explosions comes continued accretion after the explosion sets in. We show this in the bottom panel of Figure 6, where we show the PNS mass (PNS defined as the region with a matter density of at least 10¹¹ g cm⁻³) in 2D (solid lines) and 1D (dashed lines). The simulations that do not produce explosions have similar PNS masses in 1D and 2D. For the simulations that explode, the growth of the PNS mass slows after the explosion sets in, but is not zero in any model at the latest times simulated. Finally, we note that the explosions are highly nonspherical, with strong accretion occurring along various directions. We show graphical snapshots of our three exploding models in Figure 7 at a time determined by when the maximum shock radius reaches 4000 km, corresponding to 1340, 1070, and 830 ms after bounce for models s15, s20, and s25, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.3.1. Resolution and Random Perturbations

We simulate core collapse in model s25 with a lower resolution. As discussed in Section 2.3, our production simulations have a central cell size of ∼488 m and maintain an effective angular (Δr/r) resolution from at least 05 outside ∼55 km. This results in the first refinement decrement (from ∼500 m to ∼1 km) at ∼107 km. For our lower-resolution simulation, we only enforce an effective resolution of 09, giving refinement decrements at ∼65 km, ∼139 km, and so on. In addition to this simulation, we also simulate model s25 at our standard resolution and with 18 energy groups instead of our standard 12. In both tests, explosions develop at similar times to our baseline model. The development of the explosion energies for the successful models (see dashed and dot-dashed lines in Figure 6) is similar between these tests.

As a further test, we explore the effect of imposed random perturbations on the development of the explosion. We simulate model s20 five times, each time perturbing the initial density field (when we transition from our 1D model to the 2D simulations at ∼15 ms after bounce) randomly (and uniformly) by up to ±0.1%. We show the evolution of the shock radii in Figure 8. Stochastic variations in the hydrodynamic instabilities behind the shock cause the explosion times (determined by when the mean shock passes 400 km) to vary from ∼670 to ∼830 ms. The average is ∼723 ms with a formal variance of ∼63 ms. We note that the explosion time of our main simulation of s20 falls near the latter end of this range. These simulations are also discussed in O'Connor et al. (2017).

Zoom In Zoom Out Reset image size

4.1. GR versus Newtonian Gravity

From the results presented in the previous section, it is clear that the approximation of Newtonian gravity gives less favorable conditions for a successful core-collapse supernova explosion than the more realistic GR treatment. Even though we do not include full GR gravity and hydrodynamics, the test results of Section 3.1, as well as the comparisons to the previous work of Liebendörfer et al. (2005), Marek et al. (2006), and O'Connor (2015), convince us that using an effective GR potential allows us to capture the important aspects of GR for modeling core-collapse supernovae. Most importantly, the 1D results of Section 3.1 show that GR gravity leads to higher neutrino luminosities and higher rms energies, which directly follow from a more compact PNS structure. The neutrino mechanism relies on the capture of neutrinos emitted from the cooling region into the gain region. This capture is greater with higher luminosities, and more efficient with higher neutrino energies. However, in 1D, the added benefits of GR do not lead to successful explosions as has been previously demonstrated (Liebendörfer et al. 2001). This trend of increased heating and higher heating efficiencies carries over to 2D, which we demonstrate in Figure 9. In this figure, we show the net heating rate in the gain region (top), the heating efficiency (middle; Equation (15)), and the nonradial kinetic energy in the gain region (bottom) for the four WH07 models using the GR effective potential gravity (solid lines) and Newtonian gravity (dashed lines). In all models, both the heating rate and the heating efficiency are larger in the simulations using the GR effective potential, by up to ∼20%–30% between 100 and 200 ms after bounce. Furthermore, the appearance of the gain region occurs sooner for the GREP models by ∼20 ms. We also generally see more nonradial kinetic energy in the gain region, which is a sign of stronger neutrino-driven convection. As a consequence of these improved neutrino heating conditions, we see earlier explosions in the GREP models s20 and s25 by ∼0.5 s and an explosion in the GREP model s15 where the Newtonian model did not explode.

Zoom In Zoom Out Reset image size

4.2. Comparison to Literature

4.2.1. Newtonian Gravity

In Table 1 we have summarized other studies in the literature that simulate the four models simulated here and use Newtonian gravity. These include Dolence et al. (2015), Suwa et al. (2016), and Pan et al. (2016). Many aspects of our results compare quite well to the Newtonian results of Dolence et al. (2015). Dolence et al. (2015) use the same progenitor models, the same neutrino microphysics (based on Bruenn 1985; Burrows et al. 2006), and also a truly multidimensional neutrino transport scheme, albeit MGFLD. Matching our simulations, none of the 2D simulations in Dolence et al. (2015) explode before 600 ms, which is when their simulations stop. Comparison of τ_adv/τ_heat shows good agreement. At ∼200 ms, τ_adv/τ_heat ∼ 0.3–0.4 in both simulation sets and for all four models. Furthermore, our heating efficiencies (see Figure 9) near ∼200 ms after bounce are ∼0.06, similar to the values of ∼0.065 – ∼0.085 seen by Dolence et al. (2015). 2D Newtonian simulations of these models are also presented in Suwa et al. (2016) using the isotropic diffusion source approximation (IDSA) for neutrino transport. There, the LS220 EOS is used. They observe an early explosion in a few models, including s12, and failure in the majority of the rest. One issue preventing a detailed comparison is that Suwa et al. (2016) ignore heavy-lepton neutrinos in their simulations. It is likely the case that this prevents the PNS from cooling and causes the shock to initially stall at ∼240 km and hover at ∼200 km in failed models. Pan et al. (2016) also perform Newtonian axisymmetric simulations of core collapse using IDSA for neutrino transport. However, they find early explosions for all progenitors simulated: 11 M_⊙, 15 M_⊙, 21 M_⊙, and 27 M_⊙ models from Woosley et al. (2002), as well as two of the models explored here: s15 and s20. It is worth noting that later simulations in Pan et al. (2017) find no explosions in models s15 and s20 with Newtonian gravity.

4.2.2. GR Gravity

We also compare to studies in the literature that use GREP gravity. These include Summa et al. (2016) and Bruenn et al. (2016). We will first discuss the results of Summa et al. (2016). Summa et al. (2016) achieve explosions in all four models studied here. Their explosion times are ∼790, ∼620, ∼320, and ∼400 ms for models s12, s15, s20, and s25, respectively. These explosion times are earlier than ours by ∼300–500 ms. We do not see an explosion in model s12 before ∼1.5 s, when our simulation was stopped. The evolution of the mean shock radius is qualitatively similar. The initial stalling radius in our simulations is ∼150 km compared to ∼160 km in Summa et al. (2016). After this, there is a period of shock recession before the ultimate shock expansion. However, the simulations of Summa et al. (2016) have a much more pronounced shock recession. This may stem from the differences in the neutrino interactions used. Summa et al. (2016) use inelastic neutrino–nucleon scattering and absorption opacities from Burrows & Sawyer (1998, 1999), whereas we use simpler elastic rates from Bruenn (1985). They use microphysical electron capture rates on heavy nuclei for the collapse phase instead of the approximate form from Bruenn (1985) implemented here. Finally, Summa et al. (2016) include neutrino–neutrino pair scattering and annihilation (to other neutrino pairs) as additional pair processes. One important consequence of these improved opacities is the more efficient cooling of the PNS, which gives higher ν_x luminosities and faster contraction as a consequence. In O'Connor et al. (2017), we have made the first steps to improve our neutrino interactions by including a correction factor of the neutral-current scattering cross sections that mimics Burrows & Sawyer (1998, 1999). There we found earlier explosion times by ∼100–150 ms. Given these lingering differences between our simulations and those of Summa et al. (2016), a comparison with even closer underlying physics is warranted. However, we find it quite encouraging that we find relatively close agreement on many quantities, including neutrino luminosities, neutrino energies, neutrino heating rates and efficiencies, and explosion diagnostics like τ_adv/τ_heat, even with the differences between our simulation codes and simulation techniques.

In contrast, our results are in qualitative disagreement with the results of Bruenn et al. (2013, 2016; labeled here as Bruenn et al.). At 100 ms, all four models in Bruenn et al. stall at a radius of ∼200 km. However, after the shock stalls, it continues to slowly move out in radius until it ultimately begins to run away. There is no period of shock recession in any model. We note, however, that an updated 2D simulation in Lentz et al. (2015) of model s15 gives a stall radius of ∼175 km; there is still no recession of the shock front in this updated simulation. We can also directly compare neutrino quantities and heating rates with the results in Bruenn et al. The neutrino luminosities predicted from both sets of simulations agree within ∼5% (this difference is in part due to the fact that Bruenn et al. measure their quantities in the co-moving frame at 1000 km where the fluid velocities are ∼0.03c). The rms energies (right column of Figure 5) show a larger difference. Bruenn et al. predict neutrino energies up to ∼20% larger than our results at ∼100 ms after bounce. This difference cannot be explained by the difference between the co-moving frame measurements alone. It could be due to differences in the neutrino interactions (which include some of the improved rates included in Summa et al. (2016) listed above, but not the neutrino–neutrino scattering or annihilation to other neutrino pairs) or the neutrino transport (which is MGFLD, a more approximate, one-moment transport scheme using the ray-by-ray+ approximation). Roughly, the neutrino heating rates scale with both the electron-type neutrino luminosity and the electron-type mean squared neutrino energy; therefore, we expect to see the influence of the higher rms energies seen in Bruenn et al. reflected in the heating rates. Our GREP simulations have a net heating rate in the gain region at 100 ms after bounce of ∼6, ∼6, ∼10, and ∼10 B s⁻¹ for models s12, s15, s20, and s25, respectively, while the heating rates in Bruenn et al. at the same time are ∼9, ∼13, ∼20, and ∼20 B s⁻¹ for models s12, s15, s20, and s25, respectively (the heating rate at 100 ms after bounce in the updated 2D simulation of model s15 in Lentz et al. 2015 is 9 B s⁻¹). The heating rates seen in the Bruenn et al. simulations are ∼50% higher than ours. This large discrepancy at such an early post-bounce time warrants further investigation beyond comparison of literature results.

4.2.3. 1D Parameterized Models

We note that the heating efficiency of successful models is significantly less than the heating efficiency one would predict was required for an explosion in 1D. For example, O'Connor & Ott (2011) predict a critical heating efficiency of 0.172, 0.176, and 0.185 for models s15, s20, and s25, respectively. The success of these models in our simulations is not due solely to the increase in neutrino heating as a result of multidimensional instabilities; rather, the instabilities themselves are directly contributing to success of the explosion (Couch & Ott 2015). Our results, while only consisting of four models, are inconsistent with the 1D parameterized model predictions of success/failure from Ertl et al. (2016): our exploding models are predicted to fail, and our failed model is predicted to succeed. To extend this further, we have simulated four additional models, s21, s22, s23, and s24 from Woosley & Heger (2007). Based on Ertl et al. (2016), s21 is predicted to explode and the remaining models are predicted to fail. In particular, models s22 through s25 are in a region of ZAMS mass that many studies have predicted will form failed supernovae and black holes. As can be seen in Figure 10, this is not the case; similar to the pattern observed in the other models, the s21 model fails to explode, while the other three models successfully explode. We note the similarity of the shock evolution in s21 to that of s15WW95. These two models have remarkably similar presupernova structure in the inner regions, characterized by a steep gradient in the density at a compositional interface at a relatively low mass coordinate (compared to the other models explored here). The ultimate predictive power of 1D parameterized explosion models like these must be further tested (or calibrated) with multidimensional models, preferably 3D models. The impact of the multidimensional dynamics on the explosion mechanism is likely significant enough that modeling it through increased neutrino heating alone may not be sufficient.

Zoom In Zoom Out Reset image size

4.3. Impact of Proto–Neutron Star Convection

Our 2D models develop PNS convection starting at ∼60 ms after bounce. As mentioned earlier, this convection has a quantitative impact on the neutrino luminosities, particularly the heavy-lepton neutrino luminosities. We will discuss PNS convection in more detail here.

In Figure 11, we show color maps of the Brunt–Väisälä frequency, ω_BV, as a function of time and radius for the GREP model s20, for both 1D (left) and 2D (middle). While convection will not develop in our 1D model, examining ω_BV is useful nonetheless. For each time slice (angle-averaging the matter properties in the case of 2D) we compute ${\omega }_{\mathrm{BV}}^{2}$ (Foglizzo et al. 2006),

$$\begin{eqnarray}&&{\omega }_{\mathrm{BV}}^{2}=\displaystyle \frac{d\phi }{{dr}}\left(\displaystyle \frac{1}{\gamma }\displaystyle \frac{d\mathrm{ln}P}{{dr}}-\displaystyle \frac{d\mathrm{ln}\rho }{{dr}}\right),\end{eqnarray} \tag{ 18 }$$

and take

$$\begin{eqnarray}&&{\omega }_{\mathrm{BV}}=\mathrm{sign}({\omega }_{\mathrm{BV}}^{2})\sqrt{| {\omega }_{\mathrm{BV}}^{2}| },\end{eqnarray} \tag{ 19 }$$

where dϕ/dr is the radial gradient of the gravitational potential and γ is the adiabatic index. Positive values of ω_BV (blue) denote convectively stable regions, while negative values (red) denote convectively unstable regions. On each graph we also show the mean and maximum shock radius as green and magenta lines, respectively, and the average gain radius as a red line. While the Brunt–Väisälä frequency shows the susceptibility to convection, the anisotropic velocity,

$$\begin{eqnarray}&&{v}_{\mathrm{aniso}}=\sqrt{\displaystyle \frac{\langle \rho [{({v}_{r}-\langle {v}_{r}\rangle )}^{2}+{v}_{\theta }^{2}]\rangle }{\langle \rho \rangle }},\end{eqnarray} \tag{ 20 }$$

shows the actual presence of convective motions. The right panel of Figure 11 shows the anisotropic velocity in the 2D simulation of s20. Also in this panel we show the location of the heavy-lepton neutrinosphere from the 1D simulation (dashed line) and from the 2D simulation (solid line). The neutrinosphere is defined as the location where the spectrum-weighted flux factor equals 0.25. We note several important results from these graphs below.

Zoom In Zoom Out Reset image size

Before discussing PNS convection in detail, we briefly discuss the development of convection in the gain region. The 1D and 2D maps of ${\omega }_{\mathrm{BV}}^{2}$ are very similar up to ∼100 ms after bounce. Even though there is a convectively unstable region close to the shock front, due to the high advection rate of material through this region, convective seeds are unable to grow to full-scale convection before the matter is accreted out of the unstable region. By ∼100 ms the accretion rate has dropped, the shock radius has increased, and the convectively unstable region has grown so that convective motions begin to develop in earnest as seen through the rise of v_vaniso. This causes the shock radius to deviate from the spherically symmetric solution at this time. The added pressure support to the material behind the shock from these anisotropic motions is crucial to aid the neutrino mechanism and help drive an explosion (Couch & Ott 2015; Müller & Janka 2015).

Figure 11 also shows the development of a region that is unstable to PNS convection between ∼15 and ∼25 km starting as soon as ∼60 ms after bounce in the 1D model. While we just show the GREP model s20, such a region is present in all models. In the corresponding 2D simulation, the onset of PNS convection occurs at this time, as seen in the color map of v_aniso, and works to maintain ${\omega }_{\mathrm{BV}}^{2}\sim 0$. This convection has a quantitative effect on the neutrinos being emitted in our 2D simulations. PNS convection effectively transfers internal energy from deep in the PNS to slightly larger radii. This helps support the PNS and prevents the neutrinospheres from receding as they do in the corresponding 1D simulations. This predominantly impacts the heavy-lepton neutrinos, as they are emitted from the deepest regions of the star, closest to the PNS convection zone. We explicitly show the heavy-lepton neutrinosphere location in the right panel of Figure 11; the solid line is the neutrinosphere from the 2D simulations, while the dashed line is from the corresponding 1D simulation. The increased neutrinosphere radius, as well as the higher matter temperatures at the typical radius where the heavy-lepton neutrinos are emitted from, gives a larger ν_x luminosity as can be seen in Figure 5. We note that increases, which have also been attributed to PNS convection, in the ν_x luminosities are seen in 1D versus 2D models as early as Buras et al. (2006a) and more recently in Radice et al. (2017). However, it is worth mentioning that Buras et al. (2006a) see a relatively smaller impact of PNS convection on the heavy-lepton neutrino luminosity (an increase of, at most, ∼20%) than seen here or in Radice et al. (2017).

The study of core-collapse supernovae is an exceedingly multiphysics problem. If anything, decades of research into the core-collapse supernova central engine have shown that in order to create realistic simulations of these astrophysical phenomena we need to include all of the important underlying physics. This includes, but is of course not limited to, the nuclear EOS, neutrino transport and neutrino interactions, GR gravity, realistic progenitor models, resolved multidimensional (3D) magnetohydrodynamics, and more. In this paper, we have extended the capabilities of the FLASH simulation package to now include an effective treatment of GR gravity and a multispecies, multigroup, multidimensional, velocity-dependent neutrino radiation transport solver. Both of these new implementations are designed for application to the core-collapse supernova problem.

We have tested our GR implementation of gravity and our neutrino transport solver with several classic test problems, including an unstable TOV migration test (in both 1D and 2D) and the M1 shadow test. Furthermore, we perform the benchmark core-collapse test case for GR neutrino radiation transport + hydrodynamics from Liebendörfer et al. (2005), which uses a base set of neutrino interactions based on Bruenn (1985), the LS180 nuclear EOS from Lattimer & Swesty (1991), and the 15 M_⊙ progenitor from Woosley & Weaver (1995). Our results closely match those of other neutrino transport methods. We critically assess and make clear the impact of the approximations in our neutrino transport scheme with this comparison.

Extending our simulations to axisymmetry, we study core collapse in four progenitors from Woosley & Heger (2007) with both Newtonian and GR gravity for at least ∼800 ms after bounce. While the benefits that a GR treatment of gravity can have for the core-collapse supernova problem have been discussed before, ours is the first study to systematically explore the differences between Newtonian and GR gravity in multiple dimensions and across several updated progenitor models. Our results clearly and unambiguously show that GR gravity significantly aids the neutrino mechanism in reenergizing the stalled supernova shock formed in the core collapse of massive stars. We obtain explosions in three of four progenitor models from Woosley & Heger (2007) when using GR gravity. Two of the four corresponding Newtonian simulations also achieve explosions; however, they are delayed by an additional ∼500 ms. There are numerous examples in the literature of complicated multiphysics simulations of core-collapse supernova with sophisticated neutrino treatments but only Newtonian gravity. The results presented here show that the approximation of Newtonian gravity can be quite detrimental in the context of core-collapse supernova explosions and should be relaxed if accurate results are sought.

We have also showed the stochastic nature of the core-collapse problem in 2D axisymmetry by simulating one model with five different sets of random initial perturbations. For this model, we find that the explosion time can vary by ∼±100 ms. Such stochastic variations must be considered when making detailed comparisons between simulations. We have also added four additional models to our GREP simulations, s21, s22, s23, and s24 from Woosley & Heger (2007). This region is predicted from several 1D parameterized explosion studies to be susceptible to black hole formation. We do not see any evidence of this; instead, we find explosions in three of four of these models.

Our results match, in several ways, other multidimensional, energy-dependent neutrino radiation transport simulations of core-collapse supernovae, including the Newtonian gravity simulations of Dolence et al. (2015) and the GR gravity simulations of Summa et al. (2016). However, differences still exist between the results of these works and the work presented here. This could be due to differences in the neutrino and nuclear microphysics, neutrino transport methods, and/or numerical setups. We find qualitative differences between our work and the work of Bruenn et al. (2013, 2016), which, again, could be due to any of the potential differences listed above. We note specifically significant differences in the amount of the neutrino heating present during the early post-bounce phase. A close comparison of methods, input physics, and results is warranted to resolve these differences.

The authors acknowledge fruitful conversations with F. Foucart, J. A. Harris, R. Hix, H. T. Janka, E. Lentz, O. E. B. Messer, C. D. Ott, and L. F. Roberts. Computations were performed on the Zwicky cluster at Caltech, which is supported by the Sherman Fairchild Foundation and by NSF award PHY-0960291. Computations were also performed on XSEDE, which is supported by National Science Foundation grant no. ACI-154856, and the gpc supercomputer at the SciNet HPC Consortium. SciNet is funded by the Canada Foundation for Innovation under the auspices of Compute Canada, the Government of Ontario, Ontario Research Fund—Research Excellence, and the University of Toronto (Loken et al. 2010). An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under contract DE-AC02-06CH11357. Support for this work was provided by NASA through Hubble Fellowship grant No. 51344.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555.

Software: FLASH (Fryxell et al. 2000; Dubey et al. 2009), GR1D (O'Connor & Ott 2010; O'Connor 2015), Matplotlib (Hunter 2007), VisIt (Childs et al. 2012).

In Figure 12, we test our GR effective potential implementation via a TOV star migration test. We follow the initial conditions from Marek et al. (2006). Our unstable TOV stars follow a Γ = 2 polytropic EOS, P = Kρ^Γ with a value of K = 1.455 × 10⁵ [cgs]. We take an initial central density of ρ_c = 4.93 × 10¹⁵ g cm⁻³ and solve the "case A" modified TOV equations to obtain the equilibrium (but unstable) structure. We simulate 4 ms of evolution in FLASH, in both 1D (red dashed line) and 2D (blue dot-dashed). We also use the full GR code GR1D to perform the migration test (black solid line). For this full GR test, we use the standard TOV equations to construct the initial profile. In the effective potential tests, the unstable star quickly (∼0.1 ms) starts its migration away from the unstable configuration and oscillates about the stable configuration. The full GR evolution is slower, especially at the beginning, which we can understand from the time dilation effect in the full GR simulation; the value of the central lapse for this configuration begins at ∼0.27 and migrates to ∼0.65. Our results compare very well to the "case A" results and the full GR CoCoNut results of Marek et al. (2006), completing between 7 and 8 oscillation cycles in the first 4 ms of evolution for the effective potentials and ∼5 cycles in the full GR case. This convinces us that we have correctly implemented the GR effective potential in both 1D and 2D in FLASH. We have also tested stable TOV star configurations. Using a stable TOV star, with the same EOS and an initial central density of 7.5 × 10¹⁴ g cm⁻³, FLASH maintains this central density over 4 ms to better than 1% in 1D and 0.5% in 2D.

Zoom In Zoom Out Reset image size

In this appendix we present the details of the numerical techniques used to solve the neutrino transport equations (Equations (7)–(11)), the microphysics used to compute the neutrino interaction coefficients, and the coupling to the hydrodynamics.

Closure: To close the hierarchy of moment evolution equations after the first two moments, we must specify a closure relation for the Eddington tensor P^ij in terms of the two lower moments E and Fⁱ. We choose the common M1 closure and perform the calculation in the fluid rest frame in order to achieve the appropriate limiting cases. Our framework follows closely that of O'Connor (2015) and Shibata et al. (2011). In the optically thick regime the radiation in the fluid rest frame is isotropic. For this limit, we write the Eddington tensor as ${K}^{{ij}}\equiv {K}_{\mathrm{thick}}^{{ij}}=J/3{h}^{{ij}}$, where ${h}^{{ij}}=1+{W}^{2}{v}^{i}{v}^{j}$ is the projection tensor. In optically thin, free-streaming regions far from the source ${K}^{{ij}}\equiv {K}_{\mathrm{thin}}^{{ij}}=J({H}^{i}{H}^{j}/{H}^{2})$. For the total Eddington tensor in the fluid rest frame we interpolate between these two limiting regimes,

$$\begin{eqnarray}&&{K}^{{ij}}=\displaystyle \frac{3(1-\chi )}{2}{K}_{\mathrm{thick}}^{{ij}}+\displaystyle \frac{3\chi -1}{2}{K}_{\mathrm{thin}}^{{ij}}.\end{eqnarray} \tag{ 21 }$$

In these equations χ is taken (Minerbo 1978) to be

$$\begin{eqnarray}&&\chi =\displaystyle \frac{1}{3}+\displaystyle \frac{2}{15}(3{f}^{2}-{f}^{3}+3{f}^{4}),\end{eqnarray} \tag{ 22 }$$

where f ≡ (H^αH_α/J²)^1/2 is the flux factor. f is equal to 0 if the radiation is isotropic, which gives a χ = 1/3 and ${K}^{{ij}}={K}_{\mathrm{thick}}^{{ij}}$. f is 1 if the radiation is fully forward peaked in some direction; for this case, χ = 1 and ${K}^{{ij}}={K}_{\mathrm{thin}}^{{ij}}$. To get the Eddington tensor in the laboratory frame, P^ij, we take the spatial components of the neutrino stress energy tensor constructed from the fluid-frame values, ${P}^{{ij}}={T}^{{ij}}={{JW}}^{2}{v}^{i}{v}^{j}+{H}^{i}{{Wv}}^{j}+{H}^{j}{{Wv}}^{i}+{K}^{{ij}}$. In practice, we compute J and H^α via Equations (12)–(14) using a guess for P^ij constructed assuming zero velocity, and then we determine Equation (21), and finally P^ij.

Explicit Fluxes: The flux terms, both spatial and energy-space, on the left-hand side of Equations (7)–(11) are computed using the value of the radiation field at the beginning of the time step. For the energy-space fluxes we follow the methods of Müller et al. (2010). For the energy-space flux term, we use the slow-motion limit expressions from Shibata et al. (2011) and the gravitational redshift terms from Rampp & Janka (2002). Specifically, we take

$$\begin{eqnarray}&&{L}^{{ij}}=\displaystyle \frac{3\chi -1}{2}\displaystyle \frac{E\,{F}^{i}\,{F}^{j}}{{F}^{k}{F}_{k}}+\displaystyle \frac{1-\chi }{2}E{\delta }^{{ij}}\end{eqnarray} \tag{ 23 }$$

and

$$\begin{eqnarray}{N}^{{ijk}}{\partial }_{j}{v}_{k} & = & \displaystyle \frac{3\chi -1}{2}\displaystyle \frac{E\,{F}^{i}\,{F}^{j}\,{F}^{k}}{{({F}^{l}{F}_{l})}^{3/2}}{\partial }_{j}{v}_{k}\\ & & +\displaystyle \frac{3(1-\chi )}{10}({F}^{i}{\partial }_{j}{v}^{j}+{F}^{j}{\partial }_{j}{v}^{i}+{F}^{j}{\partial }_{i}{v}_{j}).\end{eqnarray} \tag{ 24 }$$

The spatial fluxes are computed using an approximate Riemann solver with corrections for the optically thick regime. We describe this in detail below. To begin, we write the spatial fluxes as a finite difference between the flux at each interface,

$$\begin{eqnarray}&&{\partial }_{x}[\alpha {r}^{m}{F}^{x}]=\displaystyle \frac{{(\alpha {r}^{m})}_{k+}{{ \mathcal F }}_{k+}-{(\alpha {r}^{m})}_{k-}{{ \mathcal F }}_{k-}}{{\rm{\Delta }}x}\end{eqnarray} \tag{ 25 }$$

and

$$\begin{eqnarray}&&{\partial }_{x}[\alpha {r}^{m}{P}^{{xi}}]=\displaystyle \frac{{(\alpha {r}^{m})}_{k+}{{ \mathcal P }}_{k+}^{i}-{(\alpha {r}^{m})}_{k-}{{ \mathcal P }}_{k-}^{i}}{{\rm{\Delta }}x},\end{eqnarray} \tag{ 26 }$$

where m is either 0, 1, or 2 and we use k± to denote the k ±1/2 interface. To obtain the interface fluxes ${{ \mathcal F }}_{k+}$ and ${{ \mathcal P }}_{k+}^{i}$, we use the standard HLLE Riemann solver for hyperbolic equations. For the flux evaluation at an interface k±, we reconstruct E and Fⁱ/E to both sides of the zone interfaces using second-order TVD (total variation diminishing) interpolation. On both sides of the interface we recompute the closure as described above to obtain the cell interface values of P_ij. The characteristic speeds for the Riemann solver are calculated in a similar way to the closure in that we interpolate between the optically thick and optically thin limits. First, for each interface (characterized here by the direction k), we determine the minimum and maximum speeds on each side of the interface in both the optically thin and optically thick limits. For the optically thick limit the choice is clear,

$$\begin{eqnarray}&&{\lambda }_{\mathrm{thick},\min }^{(k)}=-\displaystyle \frac{1}{\sqrt{3}};{\lambda }_{\mathrm{thick},\max }^{(k)}=+\displaystyle \frac{1}{\sqrt{3}}.\end{eqnarray} \tag{ 27 }$$

For the thin limit, the maximum and minimum characteristic speeds are (Shibata et al. 2011)

$$\begin{eqnarray}&&{\lambda }_{\mathrm{thin},\min /\max }^{(k)}=\min /\max \left(\displaystyle \frac{\pm {F}^{k}}{\sqrt{{F}_{i}{F}^{i}}},E\displaystyle \frac{{F}^{k}}{{F}_{i}{F}^{i}}\right).\end{eqnarray} \tag{ 28 }$$

Next, to determine the maximum and minimum speed on each side of the interface, we interpolate between the optically thick (${\lambda }_{\mathrm{thick}}^{(k)}$) and free-streaming (${\lambda }_{\mathrm{thin}}^{(k)}$) regimes via

$$\begin{eqnarray}&&{\lambda }_{\min /\max }^{(k)}=\displaystyle \frac{3(1-\chi )}{2}{\lambda }_{\mathrm{thick},\min /\max }^{(k)}+\displaystyle \frac{3\chi -1}{2}{\lambda }_{\mathrm{thin},\min /\max }^{(k)}.\end{eqnarray} \tag{ 29 }$$

The final step to determine the minimum and maximum speeds for the Riemann solver is to take ${\lambda }^{(k)+}=\max ({\lambda }_{\max }^{(k),(R)},{\lambda }_{\max }^{(k),(L)})$ and ${\lambda }^{(k)-}=\min ({\lambda }_{\min }^{(k),(R)},{\lambda }_{\min }^{(k),(L)})$, where (R) and (L) denote the right- and left-hand sides of the interface, respectively.

With the reconstructed moments and minimum and maximum characteristic speeds in hand, the HLLE Riemann solution for the fluxes at the interface is then

$$\begin{eqnarray}&&{{ \mathcal F }}_{k+,\mathrm{HLLE}}\\ &&=\,\displaystyle \frac{{\lambda }^{(k)+}{F}^{k,(L)}-{\lambda }^{(k)-}{F}^{k,(R)}+{\lambda }^{(k)+}{\lambda }^{(k)-}({E}^{(R)}-{E}^{(L)})}{{\lambda }^{(k)+}-{\lambda }^{(k)-}}\end{eqnarray} \tag{ 30 }$$

and

$$\begin{eqnarray}&&{{ \mathcal P }}_{k+,\mathrm{HLLE}}^{j}=\displaystyle \frac{{\lambda }^{(k)+}{P}_{{kj}}^{(L)}-{\lambda }^{(k)-}{P}_{{kj}}^{(R)}+{\lambda }^{(k)+}{\lambda }^{(k)-}({F}_{j}^{(R)}-{F}_{j}^{(L)})}{{\lambda }^{(k)+}-{\lambda }^{(k)-}},\end{eqnarray} \tag{ 31 }$$

where here (R) and (L) label the reconstructed (or recalculated for P^ij) moments on the right- and left-hand sides of the interface k + 1/2, respectively.

However, this solution relies on the equations being hyperbolic. This condition is violated in regions with high opacity, and numerical diffusion can dominant the HLLE Riemann solution. The transition away from hyperbolicity can be traced via the Peclet number (Pe ≡ Δx(κ_s + κ_a)). If Pe is ≳1, then we must correct the Riemann solution. We follow Audit et al. (2002) and O'Connor & Ott (2013) and ultimately take the following for the interface fluxes:

$$\begin{eqnarray}&&{{ \mathcal F }}_{k+}=a\times {{ \mathcal F }}_{k+,\mathrm{HLLE}}+(1-a)\times {{ \mathcal F }}_{\mathrm{asymptotic}}^{(k)},\end{eqnarray} \tag{ 32 }$$

where $a\equiv \tanh (1/\overline{\mathrm{Pe}})$, with $\overline{\mathrm{Pe}}$ being the geometric mean of the Peclet numbers on either side of the interface, and

$$\begin{eqnarray}&&{{ \mathcal F }}_{\mathrm{asymptotic}}^{(k)}=\displaystyle \frac{4}{3}{({v}^{k}J)}_{\mathrm{upstream}}-\displaystyle \frac{1}{3\bar{\kappa }}\displaystyle \frac{{J}^{(k+1)}-{J}^{(k)}}{{\rm{\Delta }}{x}_{k}}.\end{eqnarray} \tag{ 33 }$$

For the momentum density spatial flux we use

$$\begin{eqnarray}&&{{ \mathcal P }}_{(k+1/2)}^{j}=a\times {{ \mathcal P }}_{k+,\mathrm{HLLE}}^{j}+(1-a)\times {{ \mathcal P }}_{{\rm{k}}+,\mathrm{asymptotic}}^{j},\end{eqnarray} \tag{ 34 }$$

where

$$\begin{eqnarray}&&{{ \mathcal P }}_{{\rm{k}}+,\mathrm{asymptotic}}^{j}=({P}_{{kj}}^{(R)}+{P}_{{kj}}^{(L)})/2.\end{eqnarray} \tag{ 35 }$$

In Appendix C, we perform a 2D shadow test to ensure the correct multidimensional implementation of our transport scheme. There we show one of the beneficial features of M1, the ability to maintain neutrino momentum in the free-streaming regime. In one-moment transport schemes, like flux-limited diffusion, the flux of neutrinos is determined by the spatial gradient of the zeroth moment. In the context of core-collapse supernovae, this can cause a smoothing out of the radiation field at large distances (Ott et al. 2008).

Neutrino Interactions: The last aspect of the neutrino moment evolution equations (Equations (7)–(11)) are the neutrino interaction coefficients. These coefficients are interpolated from a table computed via NuLib (O'Connor 2015), an open-source neutrino interaction library. We include a basic set of neutrino interactions and study three independent neutrino species: electron neutrino ν_e, electron antineutrino ${\bar{\nu }}_{e}$, and a characteristic heavy-lepton neutrino ${\nu }_{x}=\{{\nu }_{\mu },{\bar{\nu }}_{\mu },{\nu }_{\tau },{\bar{\nu }}_{\tau }\}$. For electron neutrino and electron antineutrino emission processes, we include electron and positron capture on protons and neutrons, respectively, following Bruenn (1985) and Burrows et al. (2006), and include weak magnetism and nucleon recoil corrections from Horowitz (2002). We also include electron neutrino emission via electron capture on heavy nuclei via the simple formalism of Bruenn (1985). Neutrino emissivities are computed from the absorption opacities via Kirchhoff's law, which equates the emissivity to the rate of absorption of an equilibrium neutrino distribution: η = κ_aE_equil. Heavy-lepton neutrino emission from electron–positron annihilation and nucleon–nucleon bremsstrahlung is handled via an approximation that computes an effective emissivity and absorption opacity and does not require energy group coupling. This approach works well for the core-collapse problem (O'Connor 2015). We include elastic scattering on neutrons and protons, which also includes weak magnetism and recoil corrections, and elastic scattering on alpha particles and on heavy nuclei, with the latter including corrections for the heavy-ion form factor, ion–ion correlations, and electron polarization (Burrows et al. 2006).

Hydrodynamic Coupling: During each time step, we solve the neutrino moment evolution equations operator-split from the hydrodynamics. We first update the hydrodynamic variables from time step n to n + 1 using FLASH's unsplit hydrodynamics solver, the fifth-order reconstruction with WENO, and the HLLC Riemann solver. We then use the updated matter variables (ρ⁽ⁿ⁺¹⁾, T⁽ⁿ⁺¹⁾, and ${Y}_{e}^{(n+1)}$) to compute the n + 1 neutrino–matter interaction coefficients, η⁽ⁿ⁺¹⁾, ${\kappa }_{a}^{(n+1)}$, and ${\kappa }_{s}^{(n+1)}$. Our neutrino radiation evolution methods (first-order time stepping and second-order spatial reconstruction) in concert with the four guard cells needed for the hydrodynamic evolution affords us the possibility to evolve two radiation substeps for each hydrodynamic step. This greatly reduces the communication needed, which will be a limiting issue in simulations involving three dimensions. For the explicit part of each of the two radiation substeps we use the most recent radiation field variables (either the n or n + 1/2 values) to determine the fluxes (both spatial and energy-space) as described above. For each radiation substep, these fluxes are treated as an explicit source term. The remaining terms, which are all local to each grid zone and only couple the neutrino energy density and momentum density of a given group and species, are solved via implicit integration to arrive at E^(n+1/2) and ${F}^{i,(n+1/2)}$ after the first substep and E⁽ⁿ⁺¹⁾ and ${F}^{i,(n+1)}$ after the second. For these substeps we assume ${P}^{{ij},(n+1/2)}\,={p}^{{ij},(n)}{E}^{(n+1/2)}$ and ${P}^{{ij},(n+1)}={p}^{{ij},(n+1/2)}{E}^{(n+1)}$. We compute the energy, momentum, and lepton exchange with the matter for each substep and update ${v}_{i}^{(n+1)}$, ${e}^{(n+1)}$, and ${Y}_{e}^{(n+1)}$ at the end of the radiation step. For each substep, the change in ρv_i is taken as

$$\begin{eqnarray}{\rm{\Delta }}[\rho {v}_{i}] & = & -4\pi \alpha {\rm{\Delta }}t\displaystyle \sum _{\nu ,{\epsilon }_{j}}{\rm{\Delta }}\epsilon [-({\kappa }_{s}+{\kappa }_{a}){H}^{i}\\ & & +W(\eta -{\kappa }_{a}J){v}^{i}],\end{eqnarray} \tag{ 36 }$$

where the sum is over all neutrino species and neutrino energy groups. The change in the total specific energy, e_tot, is taken as

$$\begin{eqnarray}{\rm{\Delta }}[\rho {e}_{\mathrm{tot}}] & = & {\rm{\Delta }}[\rho {e}_{\mathrm{int}}]+{\rm{\Delta }}[\rho {e}_{\mathrm{kin}}]\\ & = & -4\pi \alpha {\rm{\Delta }}t\displaystyle \sum _{\nu ,{\epsilon }_{j}}{\rm{\Delta }}\epsilon [W(\eta -{\kappa }_{a}J)-({\kappa }_{a}+{\kappa }_{s}){H}^{t}].\end{eqnarray} \tag{ 37 }$$

In practice, we update v_i and e_tot and then determine e_int for use in the EOS update by subtracting off the new e_kin. We explicitly mention this because in a previous iteration of this paper (O'Connor & Couch 2015) we had mistakenly ignored the contribution to the total energy due to the changing momentum. The change in Y_e is taken as

$$\begin{eqnarray}&&{\rm{\Delta }}[\rho {Y}_{e}]=-4\pi \alpha {m}_{\mathrm{amu}}{\rm{\Delta }}t\displaystyle \sum _{\nu ,{\epsilon }_{i}}{\rm{\Delta }}\epsilon \left[{s}_{\nu }\displaystyle \frac{(\eta -{\kappa }_{a}J)}{\epsilon }\right],\end{eqnarray} \tag{ 38 }$$

where s_ν is +1 for electron neutrinos, −1 for electron antineutrinos, and 0 for heavy-lepton neutrinos. The use of an operator-split method over a coupled method between the radiation and the hydrodynamics is justified because of the small time step imposed by the explicit update of the neutrino radiation fields.

In one-moment transport schemes, like flux-limited diffusion, the flux of energy density is taken to be in the direction of the spatial gradient of the neutrino energy density. This works well in the diffusion limit, but in a free-streaming regime this assumption works to artificially source neutrino momentum and wash out potentially physical gradients. One of the advantages of an M1-like transport scheme is the ability to maintain the momentum of free-streaming neutrinos. The classic test problem to demonstrate this ability of M1 transport is the shadow test. We closely follow the 2D shadow test proposed in Just et al. (2015), with the exception that we employ axisymmetry and use cylindrical coordinates, whereas Just et al. (2015) perform their shadow test in genuine 2D Cartesian coordinates; apart from a geometric factor of r, these setups are identical. At the origin, within a spherical radius of r_s = 1.5, we have a spherical emission source with an absorption opacity of

$$\begin{eqnarray}&&{\kappa }_{a}=10\exp [-{(4r/{r}_{s})}^{2}]\end{eqnarray} \tag{ 39 }$$

and an emissivity with an arbitrary normalization. Along the radial direction, a distance of 8 from the central source, we place a purely absorbing (κ_a = 10) circle (which corresponds to a torus under axisymmetry) with a radius of 2. The grid extends to +12 in the radial direction and ±5 along the symmetry axis. The grid resolution is 0.05. For this test we do not use any mesh refinement; a side effect of this is that we have free-streaming radiation on the finest resolution level of the grid, and the first-order method for the explicit flux calculation does poorly unless we lower the CFL factor from 0.4 to 0.2, which we do for this test. We do not encounter this issue with our core-collapse simulations because of our mesh refinement, the neutrino radiation is not close to the free-streaming limit until r ≳ 100 km; by that radius we have decreased our grid resolution, giving small effective CFL factors in these outer zones.

In Figure 13, we show the radiation transport solution for the neutrino energy density. To achieve a constant value in the free-streaming regime, we scale the energy density by r². We normalize the data to range between 0 and 1. We show the solution after several light-crossing times of the domain. The shadow cast by the absorbing torus is clearly seen and demonstrates this ability of our code to maintain the neutrino momentum direction. There are slight artifacts seen as ripples at larger radii and along the coordinate axes nearer to the source that arise from the first-order explicit spatial flux integration.

Zoom In Zoom Out Reset image size