The Impact of Different Neutrino Transport Methods on Multidimensional Core-collapse Supernova Simulations
Neutrinos play a crucial role in the core-collapse supernova (CCSN) explosion mechanism. The requirement of accurately calculating the transport of neutrinos makes simulations of the CCSN mechanism extremely challenging and computationally expensive. Historically, this stiff challenge has been met by making approximations to the full transport equation. In this work, we compare CCSN simulations in one- and two-dimensions with three approximate neutrino transport schemes, each implemented in the FLASH simulation framework. We compare a two moment M1 scheme with an analytic closure (M1), the Isotropic Diffusion Source Approximation (IDSA), and the Advanced Spectral Leakage (ASL) method. We identify and discuss the advantages and disadvantages of each scheme. For each approximate transport scheme, we use identical grid setups, hydrodynamics, and gravity solvers to investigate the transport effects on supernova shock dynamics and neutrino quantities. We find that the transport scheme has a small effect on the evolution of protoneutron star (PNS) radius, PNS mass, and the mass accretion rate. The neutrino luminosities, mean energies, and shock radii have a 10-20% quantitative difference but the overall qualitative trends are fairly consistent between all three approximations. We find larger differences in the gain region properties, including the gain region mass and the net heating rate in the gain region, as well as the strength of PNS convection in the core. We investigate the progenitor, nuclear equation of state, and stochastic perturbation dependence of our simulations and find similar magnitudes of impact on key quantities. We also compare the computational expense of the various approximations.
Simulations of the core-collapse supernova (CCSN) mechanism require a challenging array of input physics, including multidimensional magnetohydrodynamics, detailed neutrino transport, involved microphysics, and general relativistic gravity. Over the roughly six-decade history of computational investigation of the CCSN mechanism, this complex mix of input physics has put this problem at the cutting-edge of computational complexity and expense. Dominating the difficulty is the requirement to accurately simulate the transport of neutrinos in this complex context. The transport of neutrinos must be followed from the diffusive, optically thick protoneutron star (PNS), through the semi-transparent region between the PNS and the stalled supernova shock, and into the completely transparent, free-streaming region beyond the shock. Solving the full seven-dimensional Boltzmann transport equation for neutrinos is almost universally too steep a challenge for high-resolution, time-dependent CCSN simulations (but see ). Approximating the full transport equation is common, though different groups working on the problem apply various different approaches.
Despite the challenge of CCSN mechanism simulations, tremendous progress has been made in recent years. After decades of research marked by cycles of promise then failure, several groups are now reporting successful neutrino-driven explosions in 2D [2, 3, 4, 5, 6, 7, 8, 9, 10]. These results often show significant quantitative, and even qualitative, differences for similar initial conditions. These various works use a variety of different hydrodynamic and neutrino transport approaches, making a direct code-to-code verification impossible. Here, we present a controlled code-to-code verification of different neutrino transport approximations commonly used in multidimensional simulations of the CCSN mechanism. We compare the two-moment explicit “M1” closure method , the isotropic diffusion source approximation (IDSA) , and the advanced spectral leakage (ASL) method , using the same simulation framework, FLASH [14, 15].
In this study, we restrict ourselves to 1D and 2D simulations. While current high-fidelity simulations in 1D only result in explosions for very low-mass iron core progenitors , it has long been understood that multidimensional effects in 2D, such as neutrino-driven convection and the standing accretion shock instability (SASI), aid shock expansion and explosion . In recent years, it has become clear that fully 3D simulations are necessary as the enforced symmetry of 2D is unduly influencing the quantitative outcomes of CCSN simulations [18, 19, 20, 21, 22].
While a controlled comparison in 3D is desirable, there is yet no detailed code comparison of CCSN simulation in even 2D (although see the recent work of [23, 24]). The careful comparison of 1D, high-fidelity CCSN simulation codes executed by  is extremely valuable, and even led to improvements in the approach for approximating general relativistic gravity now widely used . Multidimensional comparisons are challenging because of the non-linear feedback between the hydrodynamics, transport, and gravity. The development of non-linear instabilities in 2D and 3D can make it a challenge to disentangle what differences in the underlying numerical methods are leading to differences in the results. Here, we attempt to address this difficulty by carrying out 1D and 2D comparisons of different neutrino transport approximations using identical hydrodynamics, gravity solvers, EoS implementations, and computational grids. Non-linear feedback and instabilities can still magnify small differences, but using a common code for everything except neutrino transport allows us the most controlled study of the impacts of different transport methods on the overall results of CCSN simulations. This approach can shed light on the magnitude of the impact of different transport approximations in multidimensional simulations and their relative computational expense.
Our presentation of this study is as follows. In Section 2, we describe our numerical methods, including the details of the three different transport approximations we employ. In Section 3, we present our results, starting with a comparison of 1D simulations then continuing on to discuss 2D simulations. We also briefly explore the impact of different equations of state (EoS) and contrast this with the differences arising from different transport schemes. Additionally, in Section 3, we discuss the difference in overall code performance for the different approaches. We conclude in Section 4.
2.1 Hydrodynamics and gravity
We use FLASH111http://flash.uchicago.edu version 4 [14, 15] for all simulations carried out in this comparison. FLASH is a publicly available, grid-based, parallel, and multidimensional simulation framework to solve the compressible hydrodynamics equations.
The general setup is identical to what have been described in [27, 28, 8, 29, 30]. Here we summarize the key features of our numerical approach. We employ the directionally-unsplit hydrodynamics (UHD) solver with third-order piecewise parabolic method (PPM)  spatial reconstruction, a hybrid slope limiter, and a hybrid Riemann solver to calculate fluxes at zone interfaces. The hybrid Riemann solver uses the HLLC Riemann solver in smooth regions but reverts to using the more diffusive HLLE Riemann solver in zones tagged as shocks to avoid odd-even oscillations .
We use the multipole Poisson solver of Couch et al.  with a maximum multipole value for the calculation of self-gravity. In order to approximate GR effects, We replace the monopole moment of the gravitational potential with an effective GR potential based on the case A implementation described in [26, 8].
We use 1D spherical and 2D cylindrical geometries in our simulations with the PARAMESH (v.4-dev ) library for adaptive-mesh-refinement (AMR). We place the outer boundary of the domain at km and employ 9 levels of refinement. We use 5 base AMR blocks in the radial dimension for both the 1D and 2D simulations, and 10 base AMR blocks along the cylindrical -direction in 2D. Each AMR block contains 16 zones per dimension, giving a smallest zone width of 0.488 km. To save computation time, we reduce the maximum AMR level based on the distance from the center of the PNS, enforcing an effective angular resolution of . We use a radial power law profile for density and velocity as outer boundary conditions to approximate the stellar envelope rather than a pure outflow boundary condition which overestimates the mass accretion flow and can affect the shock evolution at late times .
2.2 Progenitor and nuclear equation of state
We carry out our comparison using two progenitor models from  with zero-age main sequence masses of 15 and 20 (hereafter “s15” and “s20”). The structures of these two progenitor models are quite different. The s20 progenitor model has a larger and denser silicon shell compared to the s15 model, whose density declines much faster with radius in this region. This leads initially to a larger mass accretion rate onto the PNS after bounce in s20 when compared to s15. At the silicon-oxygen interface, the s20 model has a very strong density gradient which gives a marked drop in the mass accretion rate around 200 ms after bounce. Such a sharp density drop is absent in s15, further distinguishing these models. Finally, after the accretion of the silicon-oxygen interface the density structure in s20 is such that the mass accretion rate remains fairly constant, but in s15 the mass accretion rate continues to slowly decrease.
We use the Steiner, Fischer, & Hempel (SFHo) EoS, which is tuned to fit, among other parameters, neutron star (NS) radius observations . We consider two additional EoS: the Lattimer & Swesty EoS (with incompressibility MeV, LS220) , and the Hempel & Schaffner-Bielich (HS) EoS with the DD2 parameterization, HS(DD2) . Gauging the impact of different EoS gives us a point of comparison for interpreting the magnitude of the effect of different transport approximations. Variants of the Lattimer & Swesty EoS have been widely used in the CCSN simulation community since they first became available almost 30 years ago. The LS220 EoS does not, however fulfill certain theoretical and experimental nuclear physics constraints (see, e.g., [39, 40]). The HS(DD2) EoS, on the other hand, shows good agreement with nuclear experiment about cluster formation properties , but predicts larger NS radii than SFHo. All three EoS have a maximum gravitational mass greater than 2 , as required by observations .
2.3 Neutrino transport
We compare three different neutrino transport implementations using the same hydrodynamics, gravity, and EoS, described above. We use the two-moment explicit “M1” closure method , the isotropic diffusion source approximation (IDSA) scheme , and the advanced spectral leakage (ASL) method . In this section, we briefly described each method and highlight salient points of each relevant to our comparison effort.
Our M1 transport scheme is a multidimensional, three-species, energy-dependent, approximation to Boltzmann neutrino transport. Instead of evolving the entire angle-dependent distribution function, we only evolve the first two angular moments. The zeroth angular moment represents the energy density of neutrinos within an energy bin, while the first moment represents the momentum density of neutrinos within an energy bin. We follow the formulation of [42, 11, 8]. We simulate 12 energy groups, logarithmically spaced between 1 MeV and 275 MeV. We use opacities from NuLib . Briefly, these include elastic scattering on nuclei and nucleons; charged current absorption of electron type neutrinos and anti-neutrinos on nucleons and electron type neutrinos on heavy nuclei; and thermal emission of heavy-lepton neutrinos and anti-neutrinos from electron-positron annihilation and nucleon-nucleon Bremsstrahlung. These opacities and attendant corrections, including ion-ion correlations and the heavy nucleus form factor, are based on [43, 44]. We neglect weak magnetism corrections in order to more closely match the opacity sets used in IDSA and ASL.
The neutrino moment equations are solved using standard techniques borrowed from hydrodynamics. In regions of low optical depth, the evolution equations are hyperbolic and the spatial flux between grid zones is determined using a Riemann solver. In the high optical depth limit, where the optical depth of a grid zone is greater than 1, we transition the flux determination from the Riemann solution to the asymptotic diffusion limit fluxes. To close the system of equations we assume the M1 closure for the second moment . We calculate the energy space fluxes (due to gravitational red shift and velocity gradients) explicitly. The neutrino-matter interaction source terms are treated implicitly.
A detailed description of IDSA is provided in  and  Here, we briefly review the approach reiterating the equations relevant for the present comparison. In IDSA, the distribution function of transported neutrinos is decomposed into a free-streaming component and a trapped component . These two components are evolved separately and linked by a diffusion source term . The diffusion source term is expressed as
is a non-local diffusion scalar, the emissivity, the absorptivity, the scattering opacity, and the cosine of the angle between the neutrino propagation and the radial direction. The trapped neutrino distribution is evaluated using Equation 1, 2 and the transport equation
assuming the spectral shape of the trapped component to be described by a Fermi distribution function. The diffusion scalar is solved by an explicit diffusion solver. Once is determined, the net interaction rates can be evaluated by
and hydrodynamics quantities are updated by:
where is the matter density, the baryon mass and the Planck constant. is the cooling provided by and neutrinos which is modeled by a grey leakage scheme [45, 46]. In our current IDSA solver, we only consider the transport of electron flavor neutrinos and anti-neutrinos. The streaming neutrino distribution and streaming neutrino flux for the next step can be calculated from the neutrino net interaction rates and the streaming transport equation. To couple the trapped neutrino component with matter, we introduce the trapped neutrino fraction,
The quantities and are advected with the fluid. The scaling used for the trapped neutrino energy fraction ensures the inclusion of compressional heating from trapped neutrinos via neutrino pressure. These contributions have been shown to be crucial in CCSN simulations by . Note that the current IDSA solver does not include any GR corrections to the transport equations.
The IDSA solver was first implemented in 1D coupled to the AGILE hydrodynamics code  and compared with the AGILE-BOLTZTRAN code  in the Newtonian limit. The latter solves the full Boltzmann transport equation. Good agreement of neutrino fluxes and spectra between IDSA and full Boltzmann transport was found in , but IDSA leads to a slightly larger maximum shock radius (10-20%) and a faster shock contraction.
To extend the IDSA solver to multiple dimensions, one could either solve for the diffusion scalar in multiple dimensions, but keep the streaming component isotropic , or implement the IDSA with a ‘ray-by-ray plus’ approach [49, 50]. In the latter case, the domain is decomposed in several radial directions, along which the transport problem is solved separately as in spherical symmetry, but neutrino quantities can be still advected in multiple dimensions. In this paper, we implement the IDSA solver in FLASH with the former approach, keeping the diffusion scalar in multi-D but solving the streaming component isotropically. 12 energy bins that are logarithmically spaced from 3 to 200 MeV are used in the IDSA solver. We use neutrino rates for the emission, absorption, and scattering of neutrinos off neutrons, protons and nuclei from  and nucleon-nucleon Bremsstrahlung from .
The Advanced Spectral Leakage (ASL) method [51, 13] is a three-species approximate neutrino treatment designed to model neutrinos in the context of core-collapse supernovae and compact binary mergers. It is based on previous gray leakage schemes [52, 46, 53], but in addition it carries spectral information on discretized neutrino energies and models trapped neutrino components. The spectral particle and energy emission rates are computed as a smooth interpolation between local production and diffusion rates. The former are the relevant rates in optically thin conditions, while the latter are computed based on timescale arguments and become relevant in optically thick regions. The modeling of the neutrino trapping at high densities is achieved by the solution of advection equations for and (Equations (7) and (8)) in analogy with the IDSA. The reconstruction of the neutrino distribution functions from is based on weak and thermal equilibrium arguments. Neutrino absorption rates in optically thin conditions are also included, based on the calculation of the neutrino densities of the free streaming neutrinos. ASL has been implemented in several hydrodynamics codes, both Eulerian and Lagrangian, including AGILE  in spherically symmetry, FISH  and SPHYNX  in 3D.
At the moment, ASL includes neutrino emission and absorption on free nucleons and nuclei , neutrino scattering off nucleons and nuclei in the elastic approximation , electron-positron annihilation , and nucleon-nucleon Bremsstrahlung .
In our current FLASH implementation, we evaluate the neutrino emissivities and opacities
at 12 energy bins logarithmically spaced between 3 and 300 MeV. The ASL
free-streaming component follows the ray-by-ray implementation of the gray
leakage scheme available in FLASH , where we map the
hydrodynamic grid quantities onto a set of radial rays and split the
multi-D problem into several 1D calculations. In 2D, the ray-grid consists of
radial zones, linearly spaced up to km and logarithmically spaced
up to km, and 37 uniform angular zones providing a
resolution of .
On each ray, we compute the energy-dependent optical depths and the spectral neutrino density. Locally, we map these quantities back to the hydrodynamics grid and evaluate the local source terms. The fluid energy and electron fraction, as well as the neutrino particle and energy densities, are then updated in an explicit way. The implementation presented in  is purely Newtonian. Due to the inclusion of an effective GR gravitational potential, we have developed a relativistic extension that takes into account the most relevant relativistic effects (e.g., relativistic Doppler effect and gravitational redshift) in the propagation of the free streaming neutrinos in optically thin conditions, see A. Neutrino stress in the momentum equation and from free streaming particles is not taken into account in the current implementation.
ASL contains three free parameters that require calibration: , , and . As outlined in , affects the total luminosity and the heat deposition, the neutrino energy, and the PNS cooling rate. In  these parameters are calibrated in 1D against detailed Boltzmann transport, using the AGILE-BOLTZTRAN code , a 15 zero-age main sequence mass progenitor , and the LS220 EoS. The calibrated (standard) values were , , and , where is the mass fraction of heavy nuclei. Since the implementation of ASL used in this comparison also contains GR corrections in contrast to the original implementation, we have repeated the calibration using the FLASH-M1 code in 1D, the s20 progenitor and the SFHo EoS as reference case. We have obtained , , and , comparable with the standard parameter set presented in . We have tested that differences between models employing the standard and the recalibrated parameter sets do not qualitatively change the simulation outcome for the calibration setup, but the original parameters lead to undesirable quantitative discrepancies for a detailed comparison. For instance, the maximum shock radius is reached about 10 ms later.
2.4 Initial condition
We evolve the s15 and s20 progenitors from core collapse to 15 ms post-bounce using GR1D [53, 11] and then remap the simulations to FLASH. GR1D employs the same M1 scheme as in FLASH-M1 but additionally includes inelastic process that are not included in FLASH-M1, FLASH-IDSA, and FLASH-ASL (hereafter, M1, IDSA, and ASL). This post-bounce remapping approach is similar to that employed in .
Restarting a GR1D simulation with M1 is straightforward, since M1 and GR1D share identical variables and inelastic processes are subdominant during the accretion phase after core bounce. However, in IDSA and ASL only trapped neutrinos are advected with the fluid. We, therefore, have to decouple trapped neutrinos from the total neutrinos in the initial conditions obtained from GR1D. We assume that the neutrino flux in M1 is purely from the free-streaming neutrinos and use the flux factor suggested in , where it was assumed that all neutrinos with a given energy are isotropically emitted at their last scattering neutrino sphere. Therefore,
where is the neutrino flux at energy , and is the corresponding neutrino sphere. is determined from energy-dependent opacities and is defined as the radius where the energy-dependent optical depth becomes . Once we know the distribution function of the free-streaming neutrinos, the distribution function of the trapped neutrinos is simply , and the neutrino fraction and energy can be calculated by using Equations (7) and (8).
The comparison in the following sections are done with 12 energy bins in all three transport schemes, but it should be noted that the maximum energy bin in IDSA (200 MeV) is lower than the maximum energy bin in M1 (275 MeV) and ASL (300 MeV). We have tested each transport scheme with a varying number of energy bins (from 12 up to 20 bins) in both 1D and 2D. No significant differences (other that stochastic variations for the 2D simulations) were found.
3.1 Transport comparison in spherical symmetry
We perform a series of 1D spherically symmetric simulations with M1, IDSA, and ASL. We begin with a comparison of s20 with the SFHo EoS. Figure 1 shows the time evolution of shock radius, PNS radius and mass, mass accretion rate, and the mass and net neutrino heating rate in the gain region obtained with the different transport schemes. The mass accretion rate (measured at 500 km radius) and PNS mass are nearly identical in all three transport schemes, verifying that the initial conditions, and the hydrodynamic and gravity solvers are consistent in all three schemes. M1 restarts smoothly from the GR1D and we do not observe any noticeable effects due to relaxing the model on the FLASH grid. IDSA takes about 5 ms to relax the GR1D quantities and shows some small oscillations on the shock radius and neutrino quantities. ASL takes longer ( ms) to relax the GR1D model and a more pronounced oscillation of the PNS radius (likely due to the lack of neutrino pressure) is observed (see Figure 1). Since the transition from GR1D to FLASH-M1 is smooth, we attribute these early oscillations seen in IDSA and ASL to differences in the treatment of neutrino transport.
During the first 100 ms post-bounce, M1 and IDSA show a very similar shock radius evolution that peaks at 145 km, while the ASL run has a relatively smaller shock radius evolution and peaks at 135 km. At 15 ms, the shock radius in the ASL simulation becomes comparable to the M1 simulation, and after this time the evolution is similar. On the other hand, the IDSA run gives a 15% smaller shock radius at 200 ms post-bounce, but this difference becomes smaller at late times. All three schemes give a similar PNS radius and mass, but the M1 scheme has a slightly larger PNS radius, which might be due to the different treatment of the heavy-lepton neutrinos Figure 2 shows the time evolution of neutrino luminosity and mean energy obtained with the three transport schemes. When the shock has stalled at 80 ms, the IDSA (ASL) run has the highest (lowest) electron neutrino luminosity 80 erg s (60 erg s) among the three schemes, and the M1 run lies between the IDSA and ASL. The same trend can be seen in the electron anti-neutrino luminosity. The IDSA and M1 show similar shock radius and PNS radius (which is approximately the radius of the neutrinosphere), resulting in a similar gain radius, but IDSA has a higher heating rates at early time and a lower heating rates at late time. We note that IDSA has a second bump on its electron neutrino luminosity at 200 ms. This feature does not exist in either the M1 or ASL runs, but the electron neutrino luminosity after the second bump in the IDSA matches with the M1 and ASL runs at ms. A transition between two limiting cases in the IDSA diffusion solver in Equation 1 is a possible origin of this feature.
Both IDSA and ASL show slightly higher electron neutrino and electron anti-neutrino mean energy than M1 in the first 200 ms. The difference grows to 15% after 200 ms post-bounce (see Figure 2). The neutrino mean energy in ASL is usually higher than in M1. This excess reduces to only after the drop in the accretion rate due the progenitor shell interface. This larger excess reveals the challenge for leakage schemes to model extended scattering atmospheres. We note that the leakage solver for neutrinos in IDSA does not track the mean energies and therefore these are not plotted in Figure 2.
Figure 3 shows the radial profiles of density, entropy, electron fraction, and heating rates of the three transport schemes at different post-bounce times. By 240 ms, all schemes have developed a negative entropy gradient just below the PNS radius. The simplified treatment of the diffusive regime in leakage schemes prevents an effective transport and redistribution of the heat inside the optically thick PNS. As a consequence, the dominant heavy-flavor cooling at the PNS surface produces a more pronunced entropy gradient in ASL and IDSA, compared to M1. The negative heating rates of M1 outside of the shock front at 80 ms post-bounce are due to the exchange of momentum from streaming neutrinos. It should be noticed that ASL and IDSA only include neutrino compressional heating from trapped neutrinos, but no neutrino stress from free streaming neutrinos. During the entire post bounce evolution, M1 has a larger radial extent then IDSA and ASL, which was also noted in the PNS radii panel of Figure 1. Apart from this, at 240 and 400 ms, all three schemes give very consistent radial profiles, except the IDSA run has a smaller shock radius at 240 ms.
As a final note, we compare and contrast our results to the recent work of . In the referenced comparison work, several neutrino transport codes are compared in 1D using very similar conditions to those used here, the only major difference is that in this work here we neglect weak magnetism and recoil corrections. In that work, differences between various quantities among all the codes ranged from 5% to 15%, with some excursions upwards to 50% for some select quantities at late times (e.g. neutrino energies and net heating in the gain region).
3.2 A comparison with different EoS
Recent CCSN simulations suggest that the EoS could have impact on the explodeability [34, 49], on SASI activity , on the dynamics of stellar-mass black hole formation , and on gravitational wave and neutrino signals [60, 61, 30, 62]. In order to disentangle the effects of the EoS from the neutrino transport methods, we perform 1D simulations of the s20 progenitor for all the three transport methods with two additional EoSs: LS220 and HS(DD2). The time evolution of shock radius and electron neutrino luminosity can be seen in Figure 4 for the three transport schemes we consider. The LS220 runs have a later drop in the neutrino luminosity from about 250 ms to 300 ms due to a different treatment of the low density EoS that causes a different mass accretion evolution. The simulations using the HS(DD2) EoS give the largest shock radius, followed by simulations using SFHo and LS220 EoS. The runs with SFHo (LS220) EoS have the highest (lowest) electron neutrino luminosity, respectively. The neutrino luminosity with the HS(DD2) EoS is slightly lower than that with SFHo.
All three transport schemes show the same trends while varying the EoS, suggesting that the usage of different EoS has a lower impact than the usage of different transport schemes. Therefore, the differences we discussed in the previous section do not depend on the specific choice of the EoS.
3.3 Transport comparison in cylindrical symmetry
In this section, we extend our comparison to multiple dimensions by comparing the three transport schemes via 2D cylindrically-symmetric simulations. In Figure 5, we show the same quantities for our 2D simulations as we have shown for the 1D simulations in Figure 1. The overall behavior is very similar to 1D until about 100 ms when convection begins to take hold in the gain region, breaking the spherical symmetry as visible by stronger shock oscillations, and non-zero anisotropic velocities. Up to 400 ms post bounce, none of the three models shows signs of incipient explosions.
Figure 5 reveals that both the PNS mass and accretion rate evolve similarly for all treatments since they are essentially determined by the underlying progenitor structure and gravity, neither of which is strongly impacted by the neutrino transport scheme or dimensionality. This is consistent with our results in 1D (Section 3.1). Once spherical symmetry is broken and convection becomes non-linear (after 100 ms) several of these displayed quantities begin to deviate from the 1D results. The first noticeable deviation is in the shock radius (top left panel of Figure 5), which reaches roughly the same maximum radius (135 km for ASL, 150 km for M1 and 155 km for IDSA), but then has a much slower decline. When the silicon/oxygen interface accretes through the shock at 220 ms, the shock radii are between 110-130 km, which is 30-40 km more than the value at the corresponding time in 1D. This is due to the additional dynamical pressure support and dissipation from the turbulent motions behind the shock [63, 64, 65]. In the bottom-center and bottom-right panels of Figure 5, we show the mass and the neutrino heat deposition rate in the gain region, respectively. These quantities further show the qualitative effect of multidimensional dynamics on the CCSN central engine. Compared to the analogue quantities for the 1D cases in Figure 1, we notice a slower decrease of the mass in gain region and an increased heat deposition at later times. Both of these are a result of, and also contribute to, the presence of aspherical flows in the gain region and the increased shock radius. Lastly, we note that as seen in Figure 5, the PNS radius is decreasing at a slower rate in the 2D simulations compared to the equivalent 1D simulations. This is due to the presence of PNS convection.
Figure 6 shows the evolution of neutrino luminosities and mean energies for the 2D simulations. Compared to the 1D simulations, the non-spherical accretion of turbulent material onto the PNS leads to variable signals on small timescales. This is most evident at later times, after convection has fully developed, and in the electron neutrino and anti-neutrino signals, which are emitted closer to the material in the convection zone. The heavy-lepton neutrinos originate from deeper inside in the gravitational well where the fluid motions are calmer. In 2D, convection inside the PNS increases the heat transfer from the opaque center to the surface where neutrino cooling is more efficient. This results in higher luminosities for the representative species compared to 1D models (see the top-right panel in Figure 2 and Figure 6). The M1 simulation has less of an enhancement, consistent with the milder PNS convection discussed below.
In order to understand where the differences among the neutrino treatments occur, we consider radial profiles of angular averages at different simulation times (see Figure 7). These profiles confirm that also in 2D the PNS is very similar for the different schemes. However, in contrast to the 1D case, the negative entropy and lepton gradients trigger PNS convection, which leads to a flatter entropy profile below the PNS radius. The density and entropy per nucleon (first two rows in Figure 7) compare well at small radii, where all of the matter is shocked in each transport scheme. The differences that do arise at small radii are consistent with variations of the PNS radius. The ASL simulation shows the most compact PNS, while the PNS radius in M1 and IDSA are slightly larger. The radial profiles of density and entropy can differ substantially below the shock and inside the gain region where the angular averages contain both shocked and unshocked matter at various percentages for the different schemes. A direct comparison in this regime is less straightforward.
The neutrino radiation is chiefly coupled to the matter via energy and electron fraction source terms. In the last two rows of Figure 7 we show quantities related to these source terms, i.e. the electron fraction (second last row), and the rate of energy exchange between the matter and the neutrinos (bottom row; negative values means the matter is losing energy to neutrino interactions). In general, the matter begins to deleptonize after it accretes through the shock and dissociates into free neutrons and protons. The lower shock radius for ASL at 80 ms accounts for the difference in seen there. At late times ASL tends to exhibit larger electron fractions, even greater than 0.5, close to the shock position. The in this regime is set by the deleptonization rate, but also via the neutrino heating. In the ASL simulation, the marginally but systematically larger electron neutrino luminosity enhances the rate of conversion of neutrons into protons inside the heating region. Moreover, the ray-by-ray scheme tends to enlarge the relative differences between the electron neutrino and electron anti-neutrino spectra as seen in the luminosities in Figure 6.
The neutrino energy source term reveals the location and strength of the neutrino interaction, see the last row in Figure 7. Especially during the shock expansion at 80 ms, IDSA and ASL show a very similar cooling signature below the gain radius which is located at about 100 km. As in 1D, M1 cools less inside this region. Above the gain radius up to the shock radius, M1 and ASL show a similar heating signature where IDSA deposits slightly more heat. At radii above the shock ASL and IDSA have a vanishing neutrino energy source term, but M1 also takes the neutrino pressure work on the in-falling matter into account. Furthermore, we note that comparing to 1D (Fig. 3) where M1 shows a larger radial extent in all given profiles (e.g. the rising entropy between 30–60 km at 240 ms), all schemes show a closer agreement in the radial extent of their equivalent 2D profiles (Fig. 7). Here, the strong PNS convection in ASL and IDSA lessens the difference and leads to an apparent equalization of the PNS radii among the schemes.
Figure 8 shows the evolution of the angle-averaged electron fraction , Brunt-Väisälä frequency , neutrino energy deposition, and anisotropic velocity , giving more insights about the 2D effects on the different schemes. During the early shock expansion (first 100 ms), these profiles are very similar, except a slightly noise around the neutrino sphere in the first 50 ms in ASL. When convection happens (after ms), all values close to the shock surface diverge, but the central part remains comparable for all schemes and is consistent with the previous discussion on Figure 7.
The radial profiles of Brunt-Väisälä frequency and energy deposition reveal the different heating behavior between ASL and the other schemes during the first 100 ms. It shows that for ASL there is a lag of heating between 50–100 km during the early shock expansion which explains the lower maximum shock radius. The lack of this heating feature in ASL is also confirmed by the anisotropic velocity profile. It shows that convection inside the gain region sets in 50 ms later in ASL than in the other schemes. Furthermore, the profiles of anisotropic velocity shows that IDSA and ASL tend to result in more aspherical flows inside the PNS compared to M1. This is a direct consequence of the missing energy transport in the optically thick regime typical of leakage schemes, which causes larger entropy gradients in spherically symmetric models (Section 3.1), and stronger convection in cylindrically symmetric ones. In fact, convection is even stronger in ASL compared to IDSA because in the former case leakage prescriptions are assumed for all neutrino flavors, while in the latter only for heavy flavor neutrinos. This stronger convection also leads to more noise in the energy source term inside the PNS in the ASL simulation: electron anti-neutrinos produced in trapped conditions just below the PNS radius are advected with the fluid at larger densities, where their presence is suppressed by lepton degeneracy. As a consequence, their energy is converted into matter internal energy and competes with local neutrino cooling. However, due to the high densities inside the PNS, this spurious effect does not translate into strong entropy artifacts, but rather in noise in the neutrino energy source term, as visible in Figure 6. Furthermore, the profile of anisotropic velocity reveals that ASL and IDSA evolve shock deformations very early. This is visible by the spikes in anisotropic velocity at the shock front which is a result of averages considering shocked and unshocked matter. It has also been indicated by the radial profiles in Fig. 7, where at 80 ms (first column) the profile for M1 shows a sharp discontinuity at the shock position, but the profiles of ASL and IDSA are slightly smoothed. This refers to acoustic waves which evolve from the strong PNS convection and aspherical accretion and disturb the spherical symmetry of the shock surface in these schemes.
A further confirmation of the robustness of this study was done by also varying the progenitor. Performing the same setup with s15 leads to the same overall behavior as seen in Figure 9. The s15 progenitor does not have the accretion of the shell interface at about 220 ms post bounce and the shock declines much faster which leads to a better overall agreement between all schemes.
3.4 Sensitivity study
Multidimensional simulations using spherical grids [6, 49, 50, 9] require small initial density perturbations to trigger asymmetric motion, e.g. convection. In contrast, our 2D simulations have been performed on a cylindrical grid which naturally introduces perturbations due to the spherical flow on a Cartesian grid. In order to reveal the influence of these initial perturbations on the simulation outcome, but also the specific feedback on the different transport schemes, we perform five additional simulations adding random perturbations at the level of 0.1% in the initial post-bounce conditions for each transport scheme. We adopt the 2D setup with the s20 progenitor and SFHo EoS as described in Sec. 3.3.
Figure 10 shows the spread of average shock radii for simulations with 5 different perturbation seeds for each transport scheme. It reveals that the influence of these perturbations is very low during the the shock expansion phase ( ms). The first visible deviations among the runs for a given transport scheme (colored bands) are visible for IDSA when the shock stalls ( ms). For M1 and ASL the deviations begin growing later ( ms). These deviations happened when non-spherical transient waves moving around the post-shock region. At the moment when the progenitor shell interface crosses the shock (220 ms), the sudden re-expansion of the shock further broadens the deviations. Even so much as to lead to an explosion in one of the five M1 simulations. The large deviations in the runs with ASL at the same time (220 ms) might be a ray-by-ray effect in 2D which leads to enhanced post-shock fluid motions and shock deviations as described in , but when the shock declines again, the deviations shrink and become comparable to the deviations of the non-exploding bands of the other transport schemes. For these cases, the averaged trend remains very similar to the results of the first panel in Fig. 5. As conclusion, we find that IDSA is sensitive to small perturbations at early times (100 ms), and that the ray-by-ray implementation can amplify strong asymmetric shock expansions. Inclusion of random perturbations resulted in shock revival and explosion for one M1 simulation that failed otherwise This points to the overall sensitivity of CCSN simulations to progenitor perturbations, as discussed in detail by [68, 69, 70]. This is especially true in 2D where stochastic motions can trigger shock expansion which can be very favourable for the development of an explosion.
3.5 2D code performance
In this section, we give an overview of the computational performance of the different transport schemes. Our benchmark is designed as follows: we restart the 2D setup with the s20 progenitor and SFHo EoS as described in Sec. 3.3 at 100 ms post bounce. At this time, the averaged shock radius is almost at its maximum, see fig. 5, which means that the initial AMR activity has reached an almost stable configuration. The observed range spans 100 simulation steps. We compare the ratio of core-hours spent in the neutrino treatment to the core-hours spent for solving the hydrodynamics in our simulation. The runs for the different schemes were performed on different clusters, but this ratio should robustly measure the performance of the applied scheme. Additionally from the advance in simulation time during these 100 steps, we extrapolate how many steps the simulation takes to advance one millisecond in simulation time.
Figure 11 summarizes the results. As expected, the M1 transport is the most expensive scheme per step, requiring eight times as many core-hours as a single hydrodynamic step. This is a result of evolving the first two moments of the neutrino distribution function. M1 is closely followed in expense by IDSA, which requires a factor of per hydrodynamic step. IDSA spends most of this time in solving the diffusion equation. ASL is the most approximate scheme, but has the advantage in efficiency. When running on a single node, the code spends almost the same computing time on the hydrodynamic calculation as on the neutrino scheme. Regarding the advance in simulation time, M1 and ASL show a comparable time step restriction which leads to a similar amount of steps to reach 1 ms simulation time. While the ASL time step is soley based on the CFL condition for the hydrodynamics (set by the sound speed ), M1 is set by the CFL condition for the radiation transport (set by the speed of light ). Nominally, this means a time step difference of , but since M1 performs two radiation step per hydrodynamic step, the ratio is closer to . Due to the dense regions of the PNS having a very high sound speed, this ratio is close to unity. The explicit diffusion solver in IDSA requires a much smaller time step than M1 and ASL which in the current implementation is non-adaptive leading to a constant value of 2500 steps per millisecond simulation time.
The overall 2D performance of the IDSA is worse than M1 but it should be noted that the performance of FLASH-IDSA is tuned for 3D simulations and with GPU acceleration (See B). To avoid the overhead of data copying between GPU and CPUs, we have added an additional layer of AMR block loop by doing data transfer and neutrino transport at the same time. The sequential calculation on CPUs leads to the low performance of IDSA in this particular benchmark.
With regard to the memory usage for the different neutrino transport schemes, IDSA only requires 4 additional variables on the solution per grid cell, i.e. 2 () (), which refer to the additional conservation equations for trapped neutrino particle and energy fraction, see Eqs. (7) and (8). ASL also includes the trapped component of a representative heavy flavor neutrino species , which therefore results in variables per grid cell. The trapped neutrino spectra in the IDSA and ASL are re-constructed in a AMR block level and therefore do not need to be stored in each grid cell. On the other hand, M1 carries the spectral neutrino density (scalar) and flux (vector) on the grid which leads to 4 (density + flux) 3 () (energy groups), i.e. in the case of , M1 requires 144 additional variables per grid cell. In addition to grid variables, IDSA uses variables for the spherically averaged streaming source terms, and spherically averaged thermodynamics variables to solve the streaming component, where is the number of radial zones. ASL uses 37 rays and each ray takes thermodynamics variables. These ray-by-ray variables have to be copied and synchronized for each processor but it takes much less memory than grid variables. In 2D simulations, the memory consumption of the solution usually is not a limiting factor, however it should be considered in 3D simulations where the number of grid cells significantly increases and therewith the time and memory consumption for writing checkpoints.
We have presented a series of 1D and 2D simulations of the s15 and s20 progenitors from  with the SFHo, LS220, and HS(DD2) EoS, and with three different neutrino transport schemes, including M1, IDSA, and ASL. We ran all these simulations with the publicly available code FLASH. While fixing the hydrodynamics and gravity solvers, we varied the progenitor model, nuclear EoS, and the neutrino transport scheme in order to investigate the impact of different transport methods on features of CCSN simulations.
In spherically symmetric simulations, all three transport schemes show consistent results on the evolution of the shock and neutrino quantities but with variation in certain metrics at about the 10% level. The variation we observe in 2D simulations is similar to that in 1D, but multidimensional convection leads to larger PNS radii and higher neutrino luminosities. In particular, IDSA and ASL show earlier, stronger PNS convection than M1 leading to differences in the evolution of the PNS radii and neutrino luminosities. Between transport schemes, an important difference is the prediction of neutrino luminosities and mean energies. Especially at later times, these quantities still show a large spread among the schemes. We find that convection around the PNS surface could produce an imbalance of electron and electron anti neutrinos in the ASL 2D run, giving large values of electron fraction () inside the gain region. This could be an artifact from the ray-by-ray implementation in the ASL.
When testing the sensitivity of our results to the initial progenitor profile. The differences between the transport schemes show the same trends when varying the progenitor structure and EoS. When computing resources are limited, our comparison results suggest that approximate transfer schemes can have value in their potential computational efficiency and other key factors such as nuclear EoS, turbulence, dimensionality, etc., may result in larger differences than from the neutrino transport approach. ASL runs 10 times faster than M1 and IDSA, making it possible to do a large parameter space simulations in 2D or even in 3D and giving reasonable shock dynamics. However, the ASL scheme seems to inaccurately predict the evolution, which is sensitive and important for innermost nucleosynthetic yield. In that case, one might favor the M1 scheme. The IDSA scheme lie between M1 and ASL. The IDSA runs with a slightly slower speed than M1 but is more memory efficient than M1. The memory usage could be a bottleneck for GPU programming and when moving to 3D simulations. The M1 scheme with three species neutrinos could accurately captures the changes of , and the distribution and flux of neutrinos in both opaque and transparent regions, but will cost computing time and memory.
Appendix A Relativistic corrections in FLASH-ASL
The usage of ASL in relativistic simulations or in simulations employing an effective GR gravitational potential requires the introduction of the most important relativistic corrections also in the neutrino propagation. In fact, the ASL results presented in this comparative study are qualitatively different for pure Newtonian ASL models. In the latter cases, shock revivals are observed soon after 230ms post bounce, at the occurrence of the progenitor shell interface, in 2D cylindrically symmetric models. These explosions are robust with respect to variations in the ASL free parameters and relate to systematically larger (20%) neutrino mean energies that significantly enhance the heat deposition inside the gain region. A similar effect is observed also in 1D, however the increased heating rate is not strong enough to drive an explosion in more pessimistic spherically symmetric models.
In this appendix, we present the extension of ASL we have adopted in this comparative study, which includes the gravitational redshift and the Lorentz boost between the fluid and grid reference frames. They affect radiation propagation in optically thin conditions and its absorption by the moving fluid. Since in this work we perform 1D spherically symmetric simulations and, for the 2D cylindrically symmetric models, we adopt a ray-by-ray approach for the propagation of the free streaming neutrinos, we present the relativistic extension for spherically symmetric models.
We assume a radial gauge, polar slicing metric
where is the areal radius, the lapse function relating the proper time lapse of comoving observers to the coordinate time lapse ; and the angles describing a two-sphere; , being the gravitational mass (i.e., the total energy) enclosed in a sphere of radius . The lapse function is related with the effective GR gravitational potential by:
and is obtained from the effective gravitational mass , as outlined in .
All the local quantities contained inside ASL, including neutrino source terms, are computed in the fluid reference frame (FRF), distinct from the coordinate frame (CF) associated with the metric Equation (10). The neutrino field energy in the two frames are related by a boost transformation:
where is the Lorentz factor, and the radial component of the fluid velocity as measured by an observer at constant radius. Additionally, the energy of the radiation field, climbing radially out of the gravitational well by a distance , is redshifted according to
The local spectral neutrino rates are first transformed from the FRF to the CF. To design the boost transformation at a coordinate radius for the spectral rates, we consider that the amount of emitted neutrinos (per baryonic mass) is Lorentz invariant:
where is the proper time and . The energy in the neutrino field transforms according to Equation (12):
To go from the FRF to the CF, we define and make the following ansatz about in the two frames:
The CF transformed rates are used to evolve radially the neutrino luminosities, including the gravitational redshift. This is done in an operator splitting way: first, the luminosity is evolved between two neighboring radial zones according to equation (40) in . Then the redshift correction is applied over the zone separation. We consider that, moving from to , the particle luminosity is not affected by the gravitational redshift in the CF frame:
while the energy luminosity is (cf. Equation (13)):
In analogy with the boost transformation, we define and we make the following ansatz about between and :
We solve for and by imposing Equations (17) and (18). Finally, the luminosity is locally transformed back to the FRF using the inverse of the boost transformation, Equation (16), to compute the spectral neutrino densities required to compute the local absorption rates. This procedure is applied over the entire radial profile.
Appendix B IDSA performance with GPU acceleration
Since the communication in the IDSA diffusion solver is mostly associated with neighboring zones, we have ported our IDSA solver with OpenACC for GPU acceleration. Figure 12 shows the relative computing time of the FLASH-IDSA with different dimension and block size on the Swiss supercomputer, Piz Daint. The performance is evaluated with 20 energy groups and the baseline run is using the Cray XC-30 system (with NVIDIA K20X GPU) and with 16 zones per AMR block per dimension. As discussed in Section 3.5, the speedup in 2D is worse than 3D. This is because the 2D data in an AMR block is too small to fill the GPU cores. Increasing the AMR block size from to can further improve the 2D performance (see Figure 12). In this study, we need a controlled grid step for all transport schemes to understand the transport effect. Thus, the grid setup is not tuned to the best 2D performance, and therefore GPUs are not used in the performance study in Section 3.5. The new Cray XC-50 system is about 25% faster than the original XC-30 system without using GPUs. The use of P100 GPUs on XC-50 give a speed up of 2.9 in the neutrino transport region and an overall speed up of 2.3. However even with GPU acceleration, the time step in the current IDSA solver is still restricted to 4 s due to the explicit implementation in the diffusion solver.
-  Nagakura H, Iwakami W, Furusawa S, Okawa H, Harada A, Sumiyoshi K, Yamada S, Matsufuru H and Imakura A 2017 ArXiv e-prints (Preprint 1702.01752)
-  Marek A, Janka H T and Müller E 2009 A&A 496 475–494 (Preprint 0808.4136)
-  Müller E, Janka H T and Wongwathanarat A 2012 A&A 537 A63 (Preprint 1106.6301)
-  Bruenn S W, Lentz E J, Hix W R, Mezzacappa A, Harris J A, Messer O E B, Endeve E, Blondin J M, Chertkow M A, Lingerfelt E J, Marronetti P and Yakunin K N 2016 ApJ 818 123 (Preprint 1409.5779)
-  Burrows A, Vartanyan D, Dolence J C, Skinner M A and Radice D 2016 ArXiv e-prints (Preprint 1611.05859)
-  Summa A, Hanke F, Janka H T, Melson T, Marek A and Müller B 2016 ApJ 825 6 (Preprint 1511.07871)
-  Radice D, Burrows A, Vartanyan D, Skinner M A and Dolence J C 2017 ArXiv e-prints (Preprint 1702.03927)
-  O’Connor E P and Couch S M 2018 ApJ 854 63 (Preprint 1511.07443)
-  Skinner M A, Burrows A and Dolence J C 2016 ApJ 831 81 (Preprint 1512.00113)
-  Vartanyan D, Burrows A, Radice D, Skinner M A and Dolence J 2018 ArXiv e-prints (Preprint 1801.08148)
-  O’Connor E 2015 ApJS 219 24 (Preprint 1411.7058)
-  Liebendörfer M, Whitehouse S C and Fischer T 2009 ApJ 698 1174–1190 (Preprint 0711.2929)
-  Perego A, Rosswog S, Cabezón R M, Korobkin O, Käppeli R, Arcones A and Liebendörfer M 2014 MNRAS 443 3134–3156 (Preprint 1405.6730)
-  Fryxell B, Olson K, Ricker P, Timmes F X, Zingale M, Lamb D Q, MacNeice P, Rosner R, Truran J W and Tufo H 2000 ApJS 131 273–334
-  Dubey A, Antypas K, Ganapathy M K, Reid L B, Riley K, Sheeler D, Siegel A and Weide K 2009 Parallel Computing 35 512 – 522 ISSN 0167-8191
-  Kitaura F S, Janka H T and Hillebrandt W 2006 A&A 450 345–350 (Preprint astro-ph/0512065)
-  Burrows A, Hayes J and Fryxell B A 1995 ApJ 450 830 (Preprint astro-ph/9506061)
-  Couch S M 2013 ApJ 775 35 (Preprint 1212.0010)
-  Hanke F, Müller B, Wongwathanarat A, Marek A and Janka H T 2013 ApJ 770 66 (Preprint 1303.6269)
-  Lentz E J, Bruenn S W, Hix W R, Mezzacappa A, Messer O E B, Endeve E, Blondin J M, Harris J A, Marronetti P and Yakunin K N 2015 ApJ 807 L31 (Preprint 1505.05110)
-  Melson T, Janka H T, Bollig R, Hanke F, Marek A and Müller B 2015 ApJ 808 L42 (Preprint 1504.07631)
-  Melson T, Janka H T and Marek A 2015 ApJ 801 L24 (Preprint 1501.01961)
-  Just O, Bollig R, Janka H T, Obergaulinger M, Glas R and Nagataki S 2018 ArXiv e-prints arXiv:1805.03953 (Preprint 1805.03953)
-  Cabezón R M, Pan K C, Liebendörfer M, Kuroda T, Ebinger K, Heinimann O, Thielemann F K and Perego A 2018 ArXiv e-prints (Preprint 1806.09184)
-  Liebendörfer M, Rampp M, Janka H T and Mezzacappa A 2005 ApJ 620 840–860 (Preprint astro-ph/0310662)
-  Marek A, Dimmelmeier H, Janka H T, Müller E and Buras R 2006 A&A 445 273–289 (Preprint astro-ph/0502161)
-  Couch S M, Graziani C and Flocke N 2013 ApJ 778 181 (Preprint 1307.3135)
-  Couch S M and O’Connor E P 2014 ApJ 785 123 (Preprint 1310.5728)
-  Pan K C, Liebendörfer M, Hempel M and Thielemann F K 2016 ApJ 817 72 (Preprint 1505.02513)
-  Pan K C, Liebendörfer M, Couch S M and Thielemann F K 2018 ApJ 857 13 (Preprint 1710.01690)
-  Colella P and Woodward P R 1984 Journal of Computational Physics 54 174–201
-  Quirk J J 1994 International Journal for Numerical Methods in Fluids 18 555–574
-  MacNeice P, Olson K M, Mobarry C, de Fainchtein R and Packer C 2000 Computer Physics Communications 126 330–354
-  Couch S M 2013 ApJ 765 29 (Preprint 1206.4724)
-  Woosley S E and Heger A 2007 Phys. Rep. 442 269–283 (Preprint astro-ph/0702176)
-  Steiner A W, Hempel M and Fischer T 2013 ApJ 774 17 (Preprint 1207.2184)
-  Lattimer J M and Douglas Swesty F 1991 Nuclear Physics A 535 331–376
-  Fischer T, Hempel M, Sagert I, Suwa Y and Schaffner-Bielich J 2014 European Physical Journal A 50 46 (Preprint 1307.6190)
-  Hempel M, Fischer T, Schaffner-Bielich J and Liebendörfer M 2012 ApJ 748 70 (Preprint 1108.0848)
-  Krüger T, Tews I, Hebeler K and Schwenk A 2013 Phys. Rev. C 88 025802 (Preprint 1304.2212)
-  Demorest P B, Pennucci T, Ransom S M, Roberts M S E and Hessels J W T 2010 Nature 467 1081–1083 (Preprint 1010.5788)
-  Shibata M and Taniguchi K 2011 Living Reviews in Relativity 14 6
-  Bruenn S W 1985 ApJS 58 771–841
-  Burrows A, Reddy S and Thompson T A 2006 Nuclear Physics A 777 356–394 (Preprint astro-ph/0404432)
-  Hannestad S and Raffelt G 1998 ApJ 507 339–352 (Preprint astro-ph/9711132)
-  Rosswog S and Liebendörfer M 2003 MNRAS 342 673–689 (Preprint astro-ph/0302301)
-  Mezzacappa A and Bruenn S W 1993 ApJ 410 740–760
-  Liebendörfer M, Messer O E B, Mezzacappa A, Bruenn S W, Cardall C Y and Thielemann F K 2004 ApJS 150 263–316 (Preprint astro-ph/0207036)
-  Suwa Y, Takiwaki T, Kotake K, Fischer T, Liebendörfer M and Sato K 2013 ApJ 764 99 (Preprint 1206.6101)
-  Takiwaki T, Kotake K and Suwa Y 2014 ApJ 786 83 (Preprint 1308.5755)
-  Perego A, Cabezón R M and Käppeli R 2016 ApJS 223 22 (Preprint 1511.08519)
-  Ruffert M and Janka H T 1998 A&A 338 535–555 (Preprint astro-ph/9804132)
-  O’Connor E and Ott C D 2010 Classical and Quantum Gravity 27 114103 (Preprint 0912.2393)
-  Liebendörfer M, Rosswog S and Thielemann F K 2002 ApJS 141 229–246 (Preprint astro-ph/0106539)
-  Käppeli R, Whitehouse S C, Scheidegger S, Pen U L and Liebendörfer M 2011 ApJS 195 20 (Preprint 0910.2854)
-  Cabezón R M, García-Senz D and Relaño A 2008 Journal of Computational Physics 227 8523–8540 (Preprint 0809.2755)
-  Mezzacappa A and Messer O E B 1999 Journal of Computational and Applied Mathematics 109 281–319
-  Woosley S E, Heger A and Weaver T A 2002 Reviews of Modern Physics 74 1015–1071
-  O’Connor E, Bollig R, Burrows A, Couch S, Fischer T, Janka H T, Kotake K, Lentz E J, Liebendörfer M, Messer O E B, Mezzacappa A, Takiwaki T and Vartanyan D 2018 ArXiv e-prints (Preprint 1806.04175)
-  Kuroda T, Kotake K and Takiwaki T 2016 ApJ 829 L14 (Preprint 1605.09215)
-  Richers S, Ott C D, Abdikamalov E, O’Connor E and Sullivan C 2017 Phys. Rev. D 95 063019 (Preprint 1701.02752)
-  Morozova V, Radice D, Burrows A and Vartanyan D 2018 ArXiv e-prints (Preprint 1801.01914)
-  Murphy J W, Dolence J C and Burrows A 2013 ApJ 771 52
-  Couch S M and Ott C D 2015 ApJ 799 5
-  Mabanta Q A and Murphy J W 2018 ApJ 856 22
-  Buras R, Rampp M, Janka H T and Kifonidis K 2006 A&A 447 1049–1092 (Preprint astro-ph/0507135)
-  Takiwaki T, Kotake K and Suwa Y 2012 ApJ 749 98 (Preprint 1108.3989)
-  Couch S M and Ott C D 2013 ApJ 778 L7 (Preprint 1309.2632)
-  Müller B and Janka H T 2015 MNRAS 448 2141–2174 (Preprint 1409.4783)
-  Müller B, Melson T, Heger A and Janka H T 2017 MNRAS 472 491–513 (Preprint 1705.00620)
-  Woosley S E, Almgren A, Bell J B, Glatzmaier G, Kasen D, Kerstein A R, Ma H, Nugent P, Röpke F, Sankaran V and Zingale M 2007 Type Ia supernovae Journal of Physics Conference Series (Journal of Physics Conference Series vol 78) p 012081
-  Turk M J, Smith B D, Oishi J S, Skory S, Skillman S W, Abel T and Norman M L 2011 ApJS 192 9 (Preprint 1011.3514)