© ESO 2018

1. Introduction

Core-collapse supernovae (CCSNe) are violent stellar explosions of massive stars with masses ≳8 M_⊙ at the end of their stellar evolution. Recent observations provide many details on CCSNe and their remnants, for example light curves, spectra, large-scale asymmetry, and mixing (Aschenbach et al. 1995; DeLaney et al. 2010; Milisavljevic & Fesen 2015; Boggs et al. 2015). However, there are many open questions in CCSN theory and in particular the explosion mechanism still produces explosions in numerical simulations that are sub-energetic when compared to observations (see recent reviews of Foglizzo et al. 2015; Mezzacappa et al. 2015; Janka et al. 2016; Mirizzi et al. 2016; Müller 2016; Richers et al. 2017; Burrows et al. 2018).

While it is generally believed that the neutrino-driven convection and multi-dimensional hydrodynamics instabilities are key ingredients for successful explosions, modeling CCSNe in three dimensions (3D) with Boltzmann transport is still numerically expensive (Sumiyoshi & Yamada 2012; Sumiyoshi et al. 2015; Nagakura et al. 2018) and the spatial resolution is not enough to resolve hydrodynamic instabilities and turbulence with current supercomputers (Couch & Ott 2015; Radice et al. 2016, 2018). Therefore, approximate schemes for the transport of neutrinos in multiple dimensions are still necessary, although they have to be sophisticated enough to capture the essential macro- and micro-physics.

Several two-dimensional (2D) simulations with different approximate schemes have recently been performed and investigated, however the results among different groups are less consistent than expected. For instance, the same set of simulations of the 12, 15, 20, and 25 M_⊙ progenitors shows clear differences in the shock radius evolutions during the first hundred milliseconds post-bounce between the CHIMERA code from the Oak Ridge group (Bruenn et al. 2013, 2016) and the PROMETHEUS-VERTEX code from the Garching group (Summa et al. 2016), where the former uses flux-limited diffusion for the neutrino transport and the latter uses two-moment transport with a variable Eddington factor closure. The Oak Ridge group and the Garching group use the so-called ray-by-ray-plus approach, which allocates hundreds of radial rays in the computational domain. On each ray, an independent spherically symmetric neutrino transport problem is solved, while in the optically thick regime, the neutrinos are additionally allowed to laterally advect with the fluid between the rays. The CASTRO and FORNAX simulations by the Princeton group, with flux-limited diffusion (Dolence et al. 2015) and two-moment transport (Skinner et al. 2016; Vartanyan et al. 2018) do not adopt the ray-by-ray approach and they find that the “explodability” is sensitive to the neutrino transport scheme and the detailed opacities (Burrows et al. 2018).

Liebendörfer et al. (2001); Müller et al. (2012); O’Connor & Couch (2018a) report that the higher neutrino luminosities and root mean square (rms) energies in full general relativity (GR) or effective GR calculations favor stronger supernova (SN) explosions, and O’Connor & Couch (2018a) explicitly show that the effective GR potential correction (Marek et al. 2006) could turn failed SNe into explosions in both 1D and 2D. In addition, the FLASH-M1 simulations by O’Connor & Couch (2018a) give qualitatively similar results to those of the Garching group, but ignore inelastic neutrino scattering processes. Furthermore, the Newtonian calculations from the Japanese group (Suwa et al. 2016; Nakamura et al. 2015) and the Basel group (Pan et al. 2016) with the isotropic diffusion source approximation (IDSA; Liebendörfer et al. 2009) also give different results on the shock radius evolution and explodability.

It should be noted that not only are the progenitor stars, the microphysics – such as neutrino interactions during collapse and post-bounce-, and the nuclear equations of state (EOS) key factors for the CCSN problem, but also the neutrino transport schemes, the implementation of hydrodynamics, the dimensionality of the problem, and the spatial and temporal resolution are crucial for reliable SN models. Therefore, the complexity of the problem makes it difficult to conduct a detailed comparison of SN codes among different groups, which requires computationally demanding simulations to be re-run and large amounts of 3D data to be exchanged and analyzed.

Core-collapse supernova comparisons of neutrino transport schemes were first performed for the collapse phase (Mezzacappa & Bruenn 1993a,b) with a focus on neutrino-electron scattering (NES), and were later continued to the post-bounce phase by the comparison of the VERTEX and AGILE-BOLTZTRAN codes (Liebendörfer et al. 2005). This work has recently been extended to six CCSN codes by O’Connor et al. (2018). The latter work found very consistent results, but the comparison was limited to spherically symmetric calculations. Recently, Just et al. (2018) presented a detailed 1D and 2D code comparison between the ALCAR and VERTEX codes. They discuss the influence of the ray-by-ray implementation and the impact of several different physical inputs, such as neglecting the neutrino-electron scattering, velocity, and redshift effects. However, ALCAR and VERTEX codes share similar discretization schemes and the comparison is limited to 1D and 2D. Due to the larger amount of computing time required for 3D simulations, these have only been done with a few progenitor models (Lentz et al. 2015; Kuroda et al. 2016, 2017; Melson et al. 2015a,b; Takiwaki et al. 2012, 2016; Roberts et al. 2016; Ott et al. 2017; Chan et al. 2018; O’Connor & Couch 2018b; Summa et al. 2018). A comparison of 3D models has not yet been documented. Nevertheless, the nature of CCSN is inherently 3D (DeLaney et al. 2010; Couch & Ott 2013; Boggs et al. 2015; Müller & Janka 2015; Müller et al. 2016, 2017).

In this study, we present a detailed SN code comparison in 3D. In particular, we compare the neutrino treatment with the IDSA in three different hydrodynamics codes: the uniform Cartesian grid code ELEPHANT (Käppeli et al. 2011; Liebendörfer et al. 2009), the adaptive mesh refinement (AMR) code FLASH (Fryxell et al. 2000), and the smoothed particle hydrodynamics (SPH) code SPHYNX (Cabezón et al. 2017). All three codes share the same IDSA kernel and use the same set of neutrino interactions and opacities. However, due to the different discretizations of space and time in the three hydrodynamics codes, the coupling of the IDSA to the hydrodynamic quantities had to be developed individually for each code (see Sect. 2). In addition, we compare the results of these three codes with a 3D full GR code with two-moment M1 transport: fGR1 (Kuroda et al. 2016). The 1D AGILE-BOLTZTRAN code with Boltzmann transport (Liebendörfer et al. 2004) is also used for reference calculations in spherical symmetry. Nevertheless, the differences in the discretization of space and time of the four codes in combination with the 3D degrees of freedom of the simulation data do not allow us to attribute clearly distinguishable and quantifiable effects in the results to particular code features. Rather we aim to quantify a band of uncertainty into which the results of 3D SN simulations may fall if they use different approaches for the implementation of hydrodynamics and neutrino transport. Our four methods differ most with respect to the computational domain (mesh refinement, box in a sphere, Lagrangian particles), the choice and implementation of gravity (effective GR potential for the monopole term, full GR), neutrino transport (IDSA approximation, M1 transport) and parallelization (MPI, OpenMP, or using graphics processing unit -GPU- acceleration). The codes with full neutrino transport can switch on NES, while the codes based on the IDSA cannot. The goal of this comparison is not to single out a “best code” to be used in the future (it might not be one of our four), but to further promote an approach where different codes are simultaneously used for the same physical problem. If the virtues and weaknesses of each code are well-known from previous comparisons, their results can be weighted accordingly in an uncertainty-aware interpretation of the simulated physics.

This paper is organized as follows: We describe the four SN codes and their implementations of neutrino radiation transport in the following section. In Sect. 3, we define the runs that are compared and summarize the simulation data. A detailed comparison of the most relevant physical magnitudes is presented in Sect. 4. In Sect. 5, we summarize our results and conclusions.

2. Methods and implementation

In the following section we present the four hydrodynamic codes investigated in this work. We provide an overview of their main characteristics focusing on their implementation of the Euler equations, gravity evaluation, and their coupling with the neutrino treatment.

2.1. ELEPHANT

The ELEPHANT code (ELegant and Efficient Parallel Hydrodynamics with Approximate Neutrino Transport) models the collapse and post-bounce evolution of the innermost region of a massive star in 3D. ELEPHANT is based on the 3D magneto-hydrodynamics code FISH (Käppeli et al. 2011) and implements IDSA for the neutrino transport (Liebendörfer et al. 2009; Berninger et al. 2013). The IDSA is based on two key ideas: firstly, it treats trapped neutrinos and streaming neutrinos as if they were separate particle species that are linked by a production/annihilation rate Σ. Secondly, as with flux-limited diffusion, the transition between diffusion and free streaming is handled by interpolation. The main strength of ELEPHANT is its computational simplicity and efficiency. The downside of this code is an equidistant mesh that limits the computational domain to the innermost ∼400 km of the star. The evolution of the outer layers of the progenitor star is handled by the code AGILE-IDSA (Liebendörfer et al. 2002, 2009)¹, which provides spherically symmetric data for the accreting matter at the boundary of the 3D computational domain.

With respect to the input physics, ELEPHANT only implements the basics: the Lattimer-Swesty equation of state (Lattimer & Swesty 1991) with an incompressibility parameter K = 220 MeV and the neutrino interaction rates defined in Bruenn (1985). In comparison to full Boltzmann neutrino transport in spherical symmetry (Liebendörfer et al. 2004), we found that neutrino-electron scattering of electron flavor neutrinos plays a significant role in the collapse phase, but becomes less significant during the post-bounce phase. This is consistent with the findings of Bruenn (1989), Myra & Bludman (1989), and Pan et al. (2016), however the findings of Lentz et al. (2012b) show that more modern electron capture rates on nuclei outperform neutrino-electron scattering in the collapse phase. During the collapse phase, we take the accelerated deleptonization effectively into account by using a parameterized deleptonization scheme (Liebendörfer 2005). In the post-bounce phase, we use the IDSA and neglect neutrino-electron scattering. The energy loss caused by the emission of μ/τ-neutrinos is handled by a gray leakage scheme (Sect. 2.5.2).

The hydrodynamics part of ELEPHANT solves the following set of 3D conservation equations with gravitational source terms:

(1)

(2)

(3)

The calculated unknowns, as functions of time t and space x_i, are the baryonic mass density ρ, the velocity ν_i, the electron fraction Y_e, the gravitational potential φ, and the temperature T, this latter entering the equations via the specific internal energy e(ρ,T,Y_e) and the total energy E = ρe + ρν²/2. The pressure p(ρ,T,Y_e) is provided by the equation of state. We modify the gravitational potential by general relativistic corrections following Marek et al. (2006) in order to obtain a more realistic compactness of the proto-neutron star (PNS). The factor G is the gravitational constant. As a part of the IDSA, the hydrodynamics equations additionally include equations that advect the trapped electron neutrino fractions and a multiple of the neutrino entropies, , where Z_ν represents the mean neutrino specific energy. Analogous equations are solved for the corresponding quantities of the electron antineutrinos. ELEPHANT additionally includes the operator split evolution of a magnetic field B_i = 4πb_i in the hydrodynamics part,

However, magnetic fields were ignored in the present comparison because there were no data to compare with. The neutrino transport equation is approximated by the IDSA; it separates the neutrino number density into a trapped neutrino distribution function and a streaming neutrino distribution function , and evolves these two components separately. Based on and we first construct the neutrino distribution functions assuming an equilibrium spectrum for the trapped component. Then, the diffusion equation

(4)

is solved in 3D based on the angular integral of the streaming neutrinos from the previous time step. In these equations, j_ν is the spectral neutrino emissivity, χ_ν the neutrino absorptivity and ϕ_ν includes isoenergetic scattering in the mean free path (see e.g., Bruenn 1985). The non-local diffusion term α_ν is evaluated by explicit finite differencing. All other (local) variables are iterated to convergence by an implicit Newton-Raphson scheme in each zone of the computational domain: Eq. (4) sets the net interaction rates s_ν between matter and the radiation particles, which in turn determine the updates of the electron fraction Y_e and the specific internal energy e,

(5)

(6)

(7)

The rates j_ν, χ_ν and ϕ_ν in Eq. (4) can then be adjusted to the changes in Y_e and e until convergence is achieved. In these equations m_b is the mass of a baryon, c the light speed, and h is Planck’s constant.

Assuming a stationary-state solution and neglecting observer corrections (see Liebendörfer et al. 2009 for a more detailed discussion of these assumptions) we find that the integration of the net rate over the interior of a sphere with radius R delivers a useful approximation for the neutrino number luminosity 4πR²q_ν at the surface of the sphere. Hence, for the streaming neutrino flux we obtain

(8)

In the end, we are interested in the streaming neutrino density rather than the flux. The spectral neutrino density is derived from the neutrino flux by an analytically estimated flux factor ℱ(E) (Liebendörfer et al. 2009),

This spectral neutrino density enters the diffusion Eq. (4) of the next time step. Finally, the trapped neutrino fraction Y_ν and mean energy Z_ν are updated by integrating the resulting spectral over energy. Moreover, the acceleration of matter by trapped neutrino pressure differences is applied to the velocity field,

A mathematically rigorous analysis of the IDSA with all its limiting cases is given in Berninger et al. (2013). The approximations are chosen with care and always in comparison with the more comprehensive solution of AGILE-BOLTZTRAN (Liebendörfer et al. 2004). We note that ELEPHANT includes the crucial compressional heating of the trapped neutrino gas in Eq. (2), that is an O(ν/c) observer correction in the original Boltzmann equation. Without this term, the trapped neutrino gas would not even fulfill the hydrodynamic limit.

The computational domain of ELEPHANT is decomposed into n rectangular blocks of dimension nx ⋅ ny ⋅ nz zones. Each block is assigned to one MPI-task with distributed memory and enveloped by ghost zones to support the largest stencil on the boundary zones. Moreover, the block is sliced into sheets. Depending on a flag, the operations on a sheet are either parallelized by OpenMP or sent to a GPU by OpenACC statements. The OpenACC option currently is about a factor of two faster than the OpenMP option, and is therefore the default option used in this work. Additional implementation details for ELEPHANT are given in Käppeli et al. (2011) and Liebendörfer et al. (2009).

2.2. FLASH

The IDSA solver is also implemented in the publicly available FLASH code² version 4 (Fryxell et al. 2000; Dubey et al. 2008). FLASH is a parallel, multi-dimensional hydrodynamic code based on block-structured adaptive mesh refinement (AMR). Our general setup is similar to that implemented by Couch & O’Connor (2014), but we replace the radiation transfer solver with our multi-dimensional IDSA solver that is described in the previous section and in Pan et al. (2016, 2018). In FLASH, the IDSA solver uses a similar strategy as in ELEPHANT, where we solve the diffusion source in multiple dimensions but keep the streaming component isotropic.

We use the third-order piecewise parabolic method (PPM; Colella & Woodward 1984) for spatial reconstruction, and the HLLC Riemann solver and the Hybrid slope limiter for shock-capture. The nuclear equation of state unit in FLASH incorporates the finite temperature EOS routines from O’Connor & Ott (2010) and Couch & O’Connor (2014) ³.

Our current setup supports simulations in 1D spherical coordinates, in 2D cylindrical coordinates, and in 3D Cartesian coordinates. Self-gravity is calculated by the multipole solver from Couch et al. (2013) with l_max = 16 including the same effective GR potential correction as the other codes in this work.

In addition to the GR potential correction, the main difference between the FLASH-IDSA in this paper and the IDSA implemented in Pan et al. (2016) is that we have improved the ray sampling by performing a smoother grid-ray interpolation and have improved the IDSA ray resolution from 600 zones to 1000 zones. Furthermore, in the current version of FLASH-IDSA, we disable the ray iteration in the streaming solver, but use a smaller neutrino time step (dt_ν = 4 × 10⁻⁷ s). By doing this, we obtain results that are more consistent with other SN groups (e.g., the s15 and s20 progenitor from Woosley & Heger 2007 failed to explode in 2D in Newtonian cases; Dolence et al. 2015; Summa et al. 2016; Suwa et al. 2016; Pan et al. 2017; O’Connor & Couch 2018a).

Graphics processing unit acceleration on the IDSA solver with OpenACC is also implemented in FLASH, but the parallelization strategy is different from that of ELEPHANT. Instead of sending a sheet of data, we send the whole AMR block with one layer of guard cells to a GPU all at once. While doing the calculation on GPUs, we transfer the next AMR block from host (central processing unit – CPU) to device (GPU) to avoid transfer overheads. With a block size, NXB = NYB = NZB = 16 and NE = 20, we achieve an overall speedup of 2.3 with OpenACC, where NXB, NYB, NZB are the block size in each dimension without guard cells and NE is the number of energy bins.

Analysis and visualization of FLASH simulation data are provided by the yt toolkit (Turk et al. 2011).

2.3. fGR1

The fGR1 code consists of three parts, where the evolution equations of the metric, hydrodynamics, and neutrino radiation are calculated. Each of them is solved in an operator-splitting manner, but the system evolves self-consistently as a whole, satisfying the Hamiltonian and momentum constraints. We briefly explain basic equations to be solved. We note that a geometrized unit system is used in this subsection, that is, the speed of light, the gravitational constant, and the Planck constant are set to unity: c = G = h = 1.

The space-time metric is expressed in the standard (3+1) form: ds² = −α²dt² + γ_ij(dxⁱ + βⁱdt)(dx^j + β^jdt), with α and βⁱ being the lapse and shift, respectively. The spatial metric γ_ij and its extrinsic curvature K_ij are evolved using the BSSN formulation (Shibata & Nakamura 1995; Baumgarte & Shapiro 1999) expressed in terms of the standard variables γ̃_ij, W, K, Ã_ij and Γ̃ⁱ, where γ̃_ij = ψ⁻⁴γ_ij, with ψ being the conformal factor, W(=ψ⁻²) the inverse square of the conformal factor, K = γ^ijK_ij the trace of the extrinsic curvature, Ã_ij = ψ⁻⁴(K_ij − K_γij/3), and Γ̃ⁱ = −∂_jγ̃ij (see Kuroda et al. 2012 for more details).

In the radiation hydrodynamics part, we solve the spectral neutrino transport equations, where the source terms are treated self-consistently following a procedure applying the M1 closure scheme (Shibata et al. 2011).

The total stress-energy tensor is expressed as

(9)

where and are the stress-energy tensor of fluid and energy-dependent neutrino radiation field, respectively. We note that in the equation above, summation is taken over all species of neutrinos (ν_e, ν̄_e, ν_x) with ν_x representing heavy-lepton neutrinos (i.e., ν_μ, ν_τ and their anti-particles). ε represents the neutrino energy measured in the comoving frame with the fluid. For simplicity, the neutrino flavor index is omitted below.

Introducing the radiation energy (E_(ε)), radiation flux () and radiation pressure (), measured by an Eulerian observer or (J_(ε), and ) measured in a comoving frame, can be written in covariant form as

(10)

(11)

In the above equations, n^μ = (1/α, −β^k/α) is a unit vector orthogonal to the space-like hypersurface and u^μ is the four velocity of the fluid. In the truncated moment formalism (Thorne 1981; Shibata et al. 2011), the radiation energy (E_(ε)) and radiation flux () are evolved in a conservative form and is determined by an analytic closure relation (e.g., Eq. (14)). The evolution equations for E_(ε) and are given by

(12)

and

(13)

respectively.

Here, γ is the determinant of the three metric, γ ≡ det(γ_ij), is the source term for neutrino matter interactions, and M̃^μ_(ε) is defined by M̃^μ_(ε) = M^μαβ_(ε)∇_βu_α, where denotes the third rank moment of neutrino distribution function (see Shibata et al. 2011 for the explicit expression).

By adopting the M1 closure scheme, the radiation pressure can be expressed as

(14)

where χ_(ε) represents the variable Eddington factor, and and correspond to the radiation pressure in the optically thin and thick limit, respectively. They are written in terms of J_(ε) and (Shibata et al. 2011). Following Minerbo (1978); Cernohorsky & Bludman (1994), and Obergaulinger & Janka (2011), we take the variable Eddington factor χ_(ε) as

(15)

where,

(16)

In Eq. (16), h_μν ≡ g_μν + u_μu_ν is the projection operator. By the definition of F̄_(ε) in Eq. (16), we can appropriately reproduce several important neutrino behaviors; for example, neutrino trapping in the rapidly collapsing opaque core (Kuroda et al. 2016). We simultaneously solve Eqs. (15)–(16) using the Newton method to find a converged solution of χ_(ε).

The hydrodynamic equations are written in a conservative form as

(17)

(18)

(19)

(20)

where , S_i = ρhWu_i, S_ij = ρhu_iu_j + Pγ_ij, S_k^k = γ^ijS_ij, S₀ = ρhW² − P, and ϕ = log(γ)/12. The magnitude ρ is the rest mass density, W is the Lorentz factor, h = 1 + e + P/ρ is the specific enthalpy, νⁱ = uⁱ/u^t, τ = S₀ − ρW, Y_e ≡ n_e/n_b is the electron fraction (n_X is the number density of X), e and P are the specific internal energy and pressure of matter, respectively, and m_u is the atomic mass unit. P(ρ, s, Y_e) and e(ρ, s, Y_e) are given by an EOS with s denoting the entropy per baryon. On the right-hand side of Eq. (19), Dⁱ represents the covariant derivative with respect to the three metric γ_ij.

fGR1 employs a pure MPI parallelization scheme and enforces a refluxing procedure at the refinement boundary to ensure conservation laws.

2.4. SPHYNX

The application of SPH to simulate CCSNe is not new. As early as in 1992, SPH was used to perform 2D simulations (Herant et al. 1992), highlighting the roles of the convective overturn and hydrodynamical instabilities in aiding the explosion. In 2002, SPH was used to perform the very first 3D simulations of CCSNe (Fryer & Warren 2002) and in 2004 again, this time including rotation (Fryer & Warren 2004). These works were undoubtedly remarkable landmarks and the use of SPH helped to achieve this outcome for three main reasons: firstly, SPH is intrinsi-cally devised for 3D simulations without boundaries (Lucy 1977; Gingold & Monaghan 1977); secondly, it has an adaptive spatial resolution that, in general terms, follows mass density, allowing a highly resolved PNS without presenting an unaffordable computational burden; and finally, it conserves energy and momentum (both linear and angular) by construction, which makes it very suitable for simulating rotating models, where momentum transfer is critical.

In our previous work (Perego et al. 2016) we proved the capabilities of using SPH to simulate CCSNe employing a novel spectral leakage treatment (ASL). To our knowledge, this was the first time since the early works of Herant et al. and Fryer et al. that SPH had been used to simulate this scenario, and the very first time that it was done using a spectral treatment of the neutrino component within the SPH framework. In this work, we extend the neutrino treatment, coupling the IDSA with our SPH code, SPHYNX⁴.

SPHYNX includes some of the latest upgrades that have enabled a remarkable advancement in this numerical technique during recent years. Among them are: calculating gradients using an integral approach to derivatives (IAD), which has been proven to provide more accurate derivatives (García-Senz et al. 2012; Cabezón et al. 2012; Rosswog 2015; Valdarnini 2016); using a high-order, paring-resistant family of kernels (sinc-kernels), which provides higher flexibility (Cabezón et al. 2008); and implementing novel generalized volume elements, which prevent the emergence of the tensile instability and facilitate the development of hydrodynamical instabilities (Cabezón et al. 2017).

SPHYNX solves the Euler equations, derived from a variational principle (see Rosswog 2009 and references therein). It also evaluates 3D gravity with a multipolar approximation calculated via a hierarchical tree structure that is based on the Barnes-Hut algorithm (Hernquist & Katz 1989). We have additionally included the same effective GR potential correction (Marek et al. 2006) as in ELEPHANT and FLASH to replace the monopolar term of the gravitational force. Neutrinos are handled with the IDSA treatment (Liebendörfer et al. 2009). Neutrino-electron scattering effects during collapse are taken into account by evolving the electron fraction and entropy with the parametrized deleptonization scheme of Liebendörfer (2005).

To evolve the system with SPHYNX we solve the hydrodynamical equations in Lagrangian form:

(21)

(22)

(23)

(24)

(25)

(26)

where v is the velocity vector, P_tot = P_gas + P_ν, and

(27)

The index k runs over all (anti)neutrino species. As a result, Eq. (22) includes the extra contribution of the neutrino radiation field, and therewith, the stress provided by the trapped neutrinos is directly taken into account.

From the position of the SPH particles and the equation of state, we compute the local density (ρ), the gradient of the pressure (P), the gravitational acceleration (g), and the internal energy (u). Moreover, each particle carries information regarding the electron fraction (Y_e), the neutrino abundances and energies related to the neutrino trapped component, (Y_ν, Z_ν), for all neutrino flavors. The IDSA ultimately provides the rates of change for these quantities, (Ẏ_e,I, Ẏ_ν,I, Ż_ν,I), and for the internal energy, (ė_ν,I), as denoted by the subscript I.

The abundances of electrons and neutrinos are evolved explicitly from Eqs. (25) and (26). It is worth noting that, as can be seen in Eq. (24), we implement an additional energy-like equation to evolve the energy of the trapped component of the neutrinos. In this way, we can take into account both the contribution due to the compression of the neutrino field and the rate of change provided by the IDSA, in a consistent and robust way. The details on the implementation of Eq. (24) can be found in Appendix C of Perego et al. (2016).

The IDSA is primarily a local treatment of the neutrinos. Nevertheless, it requires the spectral optical depth and the distribution function of streaming neutrinos, which are non-local quantities. Additionally, a diffusion equation for the distribution function of trapped neutrinos has to be solved. For the first two, we employ the same approach as in ELEPHANT and FLASH: we assume that the streaming component is isotropic and use spherically averaged mean free paths to integrate the optical depth. Regarding the solution of a diffusion equation, which implies second derivatives, we use the approach of Brookshaw (1985), which is less sensitive to particle distribution and has been widely used in SPH simulations with radiative transfer (see e.g., Whitehouse & Bate 2004 and Jubelgas et al. 2004 for a reformulation):

(28)

where f^t is the distribution function of the trapped component of neutrinos, λ is the mean free path, while subscripts i and j refer to a particle and its corresponding neighbors, respectively. In this way, we obtain α_i, the diffusion source for each SPH particle, as defined in Eq. (7) from Liebendörfer et al. (2009).

SPHYNX uses a hybrid MPI+OpenMP parallelization scheme. Its capabilities, implementation, and results in several tests can be found in Cabezón et al. (2017).

2.5. Implementation details

In this section we discuss the setup of the simulations in the framework of every code. We also point to the similarities and differences regarding the physics implementation, domain discretization, and spatial resolution, among the different codes.

2.5.1. Initial conditions and equation of state

In order to link to both legacy models and newer ones, we start from 1D progenitor models s15s7b2 and s15 from Woosley & Weaver (1995) and Woosley & Heger (2007), respectively. All codes begin with the same initial conditions: the onset of the collapse of a massive star with 15 M_⊙ at zero age main sequence. We map the density, temperature, electron fraction, and radial velocity profiles to the 3D computational domain. This process is adapted to each code, but the resulting initial 3D radial profiles are very similar in all codes.

In this study all codes use the Lattimer–Swesty EOS (Lattimer & Swesty 1991) with incompressibility K = 220 MeV (LS220). Although recent studies have shown that the LS220 EOS does not fulfill some theoretical and experimental nuclear constraints (Krüger et al. 2013), this is one of the most popular EOS in the SN community and has been widely used in many CCSN simulations. Moreover, the study of the effect of different EOSs in CCSN simulations is not the objective of this work.

2.5.2. Additional physics, discretization, and spatial resolution

All three Newtonian codes (ELEPHANT, FLASH, and SPHYNX) have a 3D gravity solver that has a corrected monopolar term to use the effective GR potential of Marek et al. (2006). Additionally, depending on the run, they can use the parametrized deleptonization (PD) scheme during the collapse phase to mimic the effects from neutrino-electron scattering (NES; Liebendörfer 2005). In this case, during the collapse phase, the IDSA solver is only used to update the background neutrino fraction and energies. After core bounce, the PD scheme is always turned off and the IDSA takes over. In contrast to this approximation, fGR1 calculates true NES rates and includes this process in its two-moment M1 neutrino transport scheme.

An effective GR potential correction (Case A in Marek et al. 2006) is used in all Newtonian codes to replace the monopole term in their respective gravity solvers,

(29)

where Φ_GR and Φ_NW are the GR and Newtonian gravitational potentials, respectively, Φ_eff.GR stands for the effective GR correction, and l = 0 denotes the monopole term.

The modeling of μ/τ neutrinos in the IDSA runs is handled by an inexpensive gray leakage scheme largely based on the scheme presented in the appendix of Rosswog & Liebendörfer (2003). Nevertheless, it is worth mentioning that electron neutrinos and heavy neutrinos are in fact coupled, for example by direct annihilation/creation processes, but because these cross-flavor processes are not included in our present IDSA version, electron neutrinos are decoupled from heavy neutrinos in our simulations. More generally however, caution with respect to the results of the leakage scheme is in place, especially for the longer-term evolution of the post-bounce phase. But this is not the focus of this comparison and so we feel safe to evolve them using different numerical algorithms. In our scheme, the energy loss caused by the emission of heavy flavor neutrinos is computed as smooth interpolation between diffusion and production rates. The former are valid in optically thick regions, the latter in free-streaming conditions. Neutrino elastic scattering off nucleons (Bruenn 1985) is employed as an opacity source to compute the relevant energy-averaged optical depth along the radial path. Electron-positron annihilation rates integrated over the relevant phase space are considered as the dominant production channel and are implemented according to Itoh et al. (1996). The relevant diffusion rate is calibrated against post-bounce results obtained with detailed Boltzmann neutrino transport (AGILE-BOLTZTRAN code) for a 15 M_⊙ progenitor (s15s7b2).

Regarding the energy resolution, all four codes use logarithmically spaced energy bins. The three IDSA implementations cover 3–300 MeV with 20 energy bins, while fGR1 covers 1–300 MeV with 12 energy bins. The difference in energy resolution is basically due to the more detailed neutrino transport of fGR1. Being more costly than IDSA, it partially compensates the computational load with less energy bins. Neutrino energies in fGR1 are measured in the comoving frame.

For this comparison we use a resolution of the equidistant mesh in ELEPHANT of 1 km. The 3D computational domain extends to (300 km)³ and is embedded in a larger spherically symmetric computational domain that is evolved by the spherically symmetric implicitly finite differenced code AGILE-IDSA. AGILE-IDSA provides the boundary conditions for the 3D domain at 300 km radius and is used as a limiter for the entropy at densities larger than 10¹⁴ g cm⁻³, at the very center of the PNS. The latter measure was necessary to prevent an unphysical updrift of the central entropy due to the numerical dispersion of the steep entropy gradient between the unshocked innermost part and the shock-heated outer part of the PNS in the 3D hydrodynamics code.

In FLASH, we use 3D Cartesian grids for this study. The center of the progenitor star is located at the origin of the simulation box, which includes ±5000 km in each spatial direction. The central region r ≤ 32 km has the smallest zone size of 0.488 km. Additional AMR decrements based on the distance to the origin are imposed to save computational time, giving an effective angular resolution of 0.9°−1.7°.

All SPHYNX simulations presented here use 200 000 SPH particles and cover the inner ∼1.8 M_⊙ of the progenitor star, which includes the whole iron core and part of the silicon layer. As it is characteristic for SPH codes, the spatial resolution automatically follows the density, reaching 0.477 km in the PNS. Outer layers are taken into account via an external pressure that is applied to the whole domain. The value of this pressure is taken from the 1D initial profile and it affects only the most external SPH particles.

In fGR1, the 3D computational domain is a cubic box of 15 000 km in width where nested boxes with nine refinement levels are embedded in Cartesian coordinates. To save computational time however, the initial refinement level is set to five which is increased when the maximum density reaches 5 × 10¹⁰, 10¹², 10¹³, and 3 × 10¹³ g cm⁻³ during the collapse phase. Each box contains 64³ cells and the minimum grid size near the origin is Δx = 458 m. The PNS core surface (∼10 km) and stalled shock (∼110 − 220 km) are resolved by Δx = 458 m and Δx = 7.3 km, respectively.

When rotation is included in the calculation, we implement a shellular profile, following the setup of Yokozawa et al. (2015). We discuss further details of the rotating models in Sect. 4.4.

3. Data summary, definition of acronyms of runs, and abbreviation table

CCSN simulations rely on many different aspects. For this reason, making a detailed comparison of all of them is out of the scope of this work. Nevertheless, we focus on three main pillars to compare the outcomes of our simulations: the neutrino-transport method, the gravity treatment, and the inclusion of rotation. Our runs include different combinations of these three aspects for each of the hydrodynamical codes of Sect. 2. Here we give an overview of the performed runs as well as an explanation of their names and abbreviations for ease of reference in the different comparisons.

The name of a run gives information on the aforementioned pillars and codes. To easily distinguish which code and progenitor has been used with which settings, the name of a run is twofold: the first part indicates the code and progenitor, and the second part shows the chosen setting of gravity, rotation, and (with the exception of AGILE-BOLTZTRAN) the ν-transport method. The structure of the resulting abbreviations of the performed simulations is given by

(30)

The bracket on the RHS depends on the specific setting of a run and is set as given in Table 1. Below, an example is provided to illustrate the naming generated by Eq. (30). It describes an ELEPHANT run using the s15s7b2 progenitor from Woosley & Weaver (1995) with IDSA as neutrino transport, Newtonian gravity and no rotation. The dash thereby indicates that no letter is added to the run name.

Table 1 shows an overview of abbreviations. The four 3D hydrodynamical codes (ELEPHANT, FLASH, SPHYNX and fGR1) are abbreviated by using their first capital letter. An X as code name is used as a wildcard, stating all codes with a specific configuration. We use four neutrino treatments: IDSA, IDSA with parametrized deleptonization (PD), and the M1 scheme with and without NES.

General relativity and fast rotation are only considered by their first letter in the naming code, if applied. As progenitors, we use two 15 M_⊙ models (s15s7b2 and s15) from Woosley & Weaver (1995) and Woosley & Heger (2007), respectively, referring to them by the last two digits of their publication year. Table 2 gives a complete overview of the runs performed and shall be a guide to the many run names used in this work.

4. Comparison

We present here the results of our simulations with the codes introduced in Sect. 2 and compare them to each other. In order to ease the reading of the plots, we kept the same color for the same code among all plots independently of the implemented physics or used progenitor, unless stated otherwise. Namely, ELEPHANT results are always presented in blue, FLASH in green, SPHYNX in red, and fGR1 in magenta. To show the neutrino treatment, progenitor, and gravity evaluation we follow the naming formalism introduced in Sect. 3.

4.1. Neutrino transport

The nature of the following comparison is twofold: on the one hand, we compare different approximations to neutrino transport (spherically symmetric Boltzmann solver (Liebendörfer et al. 2004) versus M1 in general relativity (Kuroda et al. 2012) versus the isotropic diffusion source approximation (IDSA; Liebendörfer et al. 2009) in FLASH, SPHYNX and ELEPHANT). On the other hand, we compare three different implementations of the same IDSA equations in the three latter codes. While we succeed to a large extent in using the same progenitor star, equation of state, and reaction rates as input physics, the implementations of the hydrodynamics equations in the three codes differ significantly. It therefore remains difficult to clearly attribute differences in the results to either the neutrino transport or the hydrodynamics part.

The differences between the results of simulations that implement the same set of equations should converge to zero if the resolution of the simulation is increased. Global astrophysical simulations however, span a dynamical range of many orders of magnitude. Thus, high resolution is excessively expensive. Often, it is neither possible nor reasonable to obtain sufficient computational resources to perform simulations that are “converged” throughout the computational domain in a mathematical sense. It is much more important to limit the deviations from an exact solution in under-resolved regimes by forcing the simulation to obey fundamental conservation laws and space- or time-averaged interaction rates, even if the spatial or temporal resolution is not sufficient throughout thus allowing a mathematically converged solution to be obtained.

Additionally, one cannot even be sure that a mathematically converged solution exists. The SN scenario involves large-scale convection of eddies between the surface of the PNS and the shock front. Small-scale turbulence connected to magnetic fields might also play a significant role. Therefore, it is likely that there always remain regions in a simulation where the evolution of small-scale features is by principle unpredictable, although it remains macroscopically constrained by conservation laws and average interaction rates. In this sense it might be in vain to strive for mathematical convergence with respect to all quantities in a simulation. Nonetheless, based on the observation of SNe one can expect that some overall physics features of the explosion remain robust with respect to perturbations: for example, the occurrence of the collapse of the stellar core, the occurrence of bounce at nuclear density, the formation of an initially rather spherical standing accretion shock, the occurrence of fluid instabilities, the expansion of the shock, the development of asymmetries, the explosion, and important features of the nucleosynthesis.

Hence, we believe that it is justifiable that astrophysics codes show slightly different results even if they are “solving” the same set of equations. They should, however, agree overall on all features that are observable and common to the majority of observed SNe.

4.1.1. Stationary state approximation in the IDSA

Matter in the gravitational collapse of a massive star has a low entropy around 1 k_B per baryon. Neutron-rich heavy nuclei dominate the scattering opacity. Due to the low ratio of neutrino energy versus the mass of a representative heavy nucleus, the scattering of trapped neutrinos is nearly elastic. A downward shift of the mean neutrino energies can still occur via inelastic neutrino-electron scattering, where both scattering partners have more similar energies. Scattering off electrons happens at a much lower rate than scattering off nuclei, but for trapped neutrinos it is efficient enough to allow more neutrinos to escape from the collapsing core due to their down-scattered energy (see e.g., Cooperstein 1988; Bruenn 1989; Bethe 1990 for a more detailed explanation). This process cannot be modeled consistently with the stationary state assumption of the IDSA. The IDSA neglects the time that neutrinos need to propagate from a location A of creation to another location B of further interaction (Liebendörfer et al. 2009). In a more recent implementation of the IDSA with time-dependent propagation of streaming neutrinos it was possible to include neutrino-electron scattering (Takiwaki et al. 2014). For similar reasons, the original implementation of the IDSA cannot adequately take into consideration the propagation of heavy neutrinos through the stellar core.

We solve these problems pragmatically, using a simple and (currently) popular method to parameterize the effects of neutrino-electron scattering during the collapse phase (Liebendörfer 2005) and a gray leakage scheme for the emission of mu-neutrinos and tau-neutrinos (see Sect. 2.5.2). Such effective treatment, once calibrated against more detailed transport schemes, can capture the most relevant cooling effect provided by the emission of heavy flavor neutrinos and the related evolution of the PNS radius with reasonable accuracy. Nevertheless, we stress that the investigation of more subtle effects related for example to the detailed treatment of neutrino-matter opacities would require a more sophisticated treatment for mu and tau neutrinos as well.

After core-bounce, the nuclei in the trapped regime become shock-dissociated. Then, the emission and absorption of neutrinos on free nucleons dominate possible effects of neutrino-electron scattering such that the approximations of the IDSA should become adequate to model the post-bounce evolution, whether neutrino-electron scattering is included or not.

The results of the collapse phase with and without parameterized deleptonization are compared in Sect. 4.2. In this section we focus on the post-bounce evolution where both transport methods, the IDSA, and full transport can be compared.

4.1.2. Comparison of data with respect to neutrino transport

In Fig. 1 the rest mass density and composition is shown at the time of bounce. The figure documents that the models start into the post-bounce phase with reasonably small initial deviations. From this point on, it is only the accretion of matter from the outer layers that controls the further evolution of the post-bounce phase. The mass accretion rate as a function of time is shown in Fig. 2 at a radius of 200 km. The accretion rate acts like a boundary condition for the evolution of the early post-bounce phase. In the Newtonian runs, the two grid-based codes show a ∼10% larger accretion than the two Lagrangian codes. However, this grouping breaks up in the general relativistic runs: SPHYNX, in GR, shows a slightly higher accretion rate than in the Newtonian case, while fGR1 seems to struggle with resolution in the early collapse phase. We note that the accretion rate evolution in SPHYNX changes slightly at 200 km when comparing Newtonian and GR-corrected simulations. This is related with the gravitational softening that is used to prevent divergences when two particles get too close. When we replace the monopolar term with the GR potential this softening is partially removed, effectively making the gravitational potential slightly stronger.

Fig. 1. Rest mass density (left column) and electron/lepton fraction profiles as a function of the enclosed baryon mass for various models at bounce. Upper and lower rows are without and with taking into account the NES effect, respectively.

Fig. 2. Time evolution of mass accretion rates measured at R = 200 km. Different colors represent simulations with different codes. Each column shows a set of simulations (see Table 2).

A broad overview of the time evolution of several key quantities during the post-bounce phase is shown in Fig. 3. The first row shows Newtonian runs with reduced input physics (labeled I in Table 1). The second row shows Newtonian runs with the effect of neutrino-electron scattering (labeled S or P in Table 1). The third and fourth rows show the corresponding general relativistic runs (labeled IG or MG with reduced physics and SG, PG or MSG with neutrino-electron scattering). The panels in the left column of the figure show the shock position as a function of time. All four panels demonstrate that within 10 ms the shock expands from a radius of a few kilometers to one of about 100 km.

Fig. 3. Time evolution of averaged shock radius, PNS radius, neutrino luminosities, and mean energies for different runs. Each row describes a comparison of different physics input as described in Table 2. Different color represents a simulation with a different hydrodynamics code.

However, the four panels in the left column also show clear differences between the runs. We believe that they are consistent with the following arguments: In the top panel it can be seen that run B95 reaches a larger shock radius than the other codes, as well as in the third panel from above, B95-G and G95-MG reach a larger shock radius than the other codes. This can be expected because in the IDSA, neutrinos are assumed to propagate instantly through the stellar core due to the steady-state approximation for streaming neutrinos. In the true transport codes B95-G and G95-MG this is not the case. There, the neutrinos require a few milliseconds to move away from the location where they are produced. During these 5 − 10 ms they block the phase space for further neutrinos to be created. Hence, the shock in runs B95-G and G95-MG suffers less from instantaneous neutrino loss than the shock in the IDSA runs. This explanation is consistent with the fact that B95 also reaches the largest shock position in the topmost panel (there is no corresponding G95-M run). It is also consistent with the observation that this difference in the early evolution of the shock radius is not seen in the second and fourth panel from the top. These panels show graphs for runs that implement the effects of neutrino-electron scattering. In these runs, the neutrino energy is lowered by down-scattering. Thus, the neutrinos are emitted earlier. In this case, the codes produce a weaker shock from the beginning and a relatively consistent shock expansion up to 10 ms for all runs.

Further deviations appear at ∼15 ms post-bounce. ELEPHANT produces a broad hump in the top panel, which can also be seen in the FLASH data, but to a much smaller extent. In the second panel from the top, both ELEPHANT and FLASH show a similar hump at a slightly earlier time. This is also present in the third panel, albeit again smaller, but cannot be distinguished in the fourth panel. We explain this hump with an early optimistic expansion of matter supported by grid alignment effects. After this first matter overturn, the alignment benefit vanishes and the shock falls back to a more steady position. In earlier, slowly rotating models these humps are less visible (compare also to runs IR and PR in Fig. 8). In SPHYNX runs, the shock evolution shows smaller humps or no humps at all. As SPH codes have no underlying structured mesh, alignment effects are suppressed. Shock oscillations, such as those present in runs S95-I and S95-P, can be linked to early prompt convection. After this first overshoot, the shocks in all models reach consistent positions with deviations below 10% between B95, E95, F95 and S95.

The second column of the panels in Fig. 3 presents the luminosity as a function of time for the electron flavor neutrinos (solid lines) and antineutrinos (dashed lines). After the peak luminosity in the neutrino burst has occurred at ∼5 ms after bounce, all luminosities decay to a value between 4 × 10⁵² erg and 6 × 10⁵² erg at 50 ms post-bounce. At first glance, the neutrino luminosities seem to oscillate without correlation during this decay. Closer inspection reveals that there is a strong correlation between the dynamics of the shock front discussed in the previous paragraph and the evolution of the electron neutrino luminosity. For example, run E95-I shows a prominent recession of the shock at 24 ms post-bounce. At exactly this time the corresponding neutrino luminosity exhibits an equally prominent peak that is emitted by electron capture when the hot matter behind the shock is compressed as it slams into the surface of the PNS. The same feature appears in runs E95-P and S95-P at 16 ms post-bounce, in run S95-I at 20 ms post-bounce, in run F95-P at 18 ms post-bounce, and in run G95-MG at 24 ms post-bounce. Most features in the neutrino luminosity can be explained in this way by the different dynamics of the shock.

One exception is perhaps the luminosity of run B95: the prominent recession of the shock at 12 ms post-bounce would suggest a rather high neutrino luminosity at that time. The luminosity of B95 does show a local maximum, but it does not exceed the luminosity of E95-I and F95-I. Here one should additionally consider the PNS radius (dashed lines in the leftmost panels), which initially shows a small PNS radius and a large shock radius. In this most optimistic run, a substantial mass hovers between these radii and has not yet settled on the PNS. At 26 ms post-bounce, the situation is reversed: as soon as the shock has receded after 14 ms post-bounce the PNS radius becomes even larger than in the other runs. At this time, the luminosity of E95-I consistently exceeds the luminosities of the other runs. The luminosity of run S95-I remains mostly below the luminosities of the other runs because the shock expands slowly and exhibits no recession except once around 20 ms post-bounce.

The antineutrino luminosities depend much less on the accretion rate. We find agreement of about 10% between the runs that use the IDSA and agreement after ∼20 ms post-bounce between the two runs that implement more sophisticated neutrino transport. However, the IDSA-runs show consistently higher antineutrino luminosities than the transport-runs. We note that a similar deviation of the antineutrinos in the IDSA has been noticed independently in O’Connor et al. (2018), but not in Liebendörfer et al. (2009). We also do not yet know why the antineutrino luminosities calculated by the two transport codes B95 and fGR1 differ by about a factor of two between 10 ms and 20 ms post-bounce. These questions certainly merit further investigation and demonstrate how important such comparisons are.

The mean energy of the neutrino flux, calculated by dividing the energy luminosity by the number luminosity, is shown in the rightmost column of the panels. In runs X95-I and X95-IG we find satisfactory agreement with deviations of the order of 20% during the most dynamic phase and of the order of 10% after the shock has assumed a more steady expansion. The longer-term neutrino and antineutrino energies in ELEPHANT are consistently higher than in FLASH and SPHYNX while the neutrino and antineutrino luminosities in fGR1 are significantly lower than in the other codes. The neutrino energies of ELEPHANT agree well with those of AGILE-BOLTZTRAN in the runs without neutrino-electron scattering and are higher than those of AGILE-BOLTZTRAN in the runs that include neutrino-electron scattering. This is consistent with the fact that with the IDSA the (anti-)neutrinos leave the star with unchanged mean energy, while they are able to downscatter in the runs with AGILE-BOLTZTRAN and fGR1, which both implement full neutrino-electron scattering.

We discover the most prominent differences between the codes in the luminosities of the μ- and τ-neutrinos. The energy loss due to the emission of heavy neutrinos in the Newtonian runs, which are based on a gray leakage scheme for the heavy neutrinos, is smaller than the corresponding energy loss calculated in AGILE-BOLTZTRAN with full neutrino transport by almost a factor two. Because the AGILE-BOLTZTRAN results assume similar values to those of an earlier comparison (Liebendörfer 2005), we assume that the AGILE-BOLTZTRAN data are correct and that the gray leakage scheme is too simple for the conditions that prevail in Newtonian runs. This is consistent with the observation that the results of fGR1 are closer to the results of AGILE-BOLTZTRAN than to the results of the leakage scheme. In the more compact PNS of the general relativistic runs, where higher temperatures are reached, the leakage scheme performs better and approaches the AGILE-BOLTZTRAN results to around 20% accuracy.

The differences in the energy loss rate by the emission of heavy neutrinos points to the possibility that some differences in Fig. 3 could be caused by a different evolution of the PNS. Figure 4 shows the central density as a function of time. At first we note that the central density of AGILE-BOLTZTRAN (orange dashed line) appears to have relatively low values. This is a result of the adaptive mesh, which pulls grid points from the center of the PNS away to the shock front. The center of the PNS is then less resolved than in the other codes. As the innermost value in AGILE-BOLTZTRAN is a cell-centered cell average, it is lower than the actual central density at the inner edge of the cell. In order to cure this inconsistency in the comparison, we extrapolate the density profile of AGILE-BOLTZTRAN to the center of the computational domain using the constraints that the density profile around the center be quadratic with zero derivative at r = 0. Furthermore, it must exactly reproduce the density averages in the first and second cells given by the output files of AGILE-BOLTZTRAN. After the extrapolation to the center, the density of AGILE-BOLTZTRAN (orange solid line) shows better agreement with the central density of the other codes.

Fig. 4. Central density evolution in all codes within a Newtonian (left panel) and GR (right panel) framework. The dashed line shows AGILE-BOLTZTRAN data before the extrapolation to the center (cf. text).

In the Newtonian runs (left panel) we find an agreement to better than 5% among all codes. However, the agreement in the general relativistic runs is less perfect. E95-PG produces the most compact PNS with the largest central density. G95-MSG produces the PNS with the lowest central density. The central density depends on the self-gravity of the neutron star, the central entropy and the central electron fraction. We find good agreement with respect to the electron fraction. The Y_e is set by the parameterized deleptonization scheme and evolves in all models in a similar fashion. Differences are found in the central entropy. Consistent with the density differences, E95-PG shows the lowest entropy of 1.2 k_B per baryon while G95-MSG produces 1.4 − 1.5 k_B per baryon. Figure 5 shows the evolution of the central entropy as a function of density during collapse (left panel) and as a function of time (right panel). The left panel shows that the main entropy difference is induced during the collapse phase when the density crosses ∼10¹² g cm⁻³. The right panel shows an unphysical slight entropy increase after bounce in the multidimensional codes that is most likely due to advection of entropy from the shock-heated layers of the PNS towards the cold center. This effect is not visible in the data of AGILE-BOLTZTRAN because the core in spherical symmetry does not require substantial advection. It is also not visible in the data of SPHYNX because advection is not an issue in the inherently Lagrangian code. The effect is smaller in ELEPHANT than in the two other 3D grid-based codes simply because the central entropy in ELEPHANT is limited by the central entropy of AGILE-IDSA to mitigate this undesired central entropy increase.

Fig. 5. Central entropy as a function of central density (left panel) and as a function of time (right panel) for GR models including NES.

If we now go back to Fig. 3 with these central entropy differences in mind, we confirm in the lower left panel that the PNS radius of G95-MSG is indeed larger than the one of E95-PG. Also the shock radius of G95-MSG is larger than the one of E95-PG. As expected, we find the inverse ordering in the neutrino luminosities: the electron flavor luminosities of G95-MSG are significantly lower than those of E95-PG. The difference of the mean energies between these two runs can be explained by the same difference of the PNS structure: the more compact PNS of E95-PG has neutrinospheres deeper in the gravitational well and produces higher mean energies than the less compact PNS of G95-MSG.

In summary, we find a consistent overall evolution in all runs. Differences during the collapse phase lead to slight differences in the entropy and compactness of the PNS. The implementations of hydrodynamics in our four codes optimize different aspects of the SN dynamics: avoidance of grid-effects in SPH, high resolution and efficiency at the shock front, adaptive capabilities and improved resolution of the PNS, and last but not least, the possibility to include full GR. These differences lead to deviations in the details of the early shock expansion. Many differences in the evolution of the neutrino luminosities can directly be linked to these hydrodynamic differences. Some deviations, nota bene, require further investigation before we can extend this comparison to later post-bounce times: Why are the electron antineutrino luminosities produced by the IDSA rather large? And under which PNS conditions is the leakage scheme for the μ- and τ-neutrinos sufficient – or insufficient?

4.2. Neutrino electron scattering

Inelastic scattering of neutrinos off electrons plays an important role in helping neutrinos to escape more freely from the center of the star as a consequence of down-scattering; it enhances deleptonization of the central core and thus leads to lower electron fraction Y_e. The electron fraction is the most important quantity to determine the stability and the size of the inner core during stellar collapse. The latter is crucial, since it defines the location of the shock at core bounce and affects the subsequent evolution of the neutrinospheres, which may be associated with the neutrino heating efficiency. The deleptonization process is mainly controlled by electron captures on free protons and elastic/inelastic scattering of electron neutrinos. In this section, we compare its treatment used in each code and discuss their properties.

In Fig. 6, we plot the central evolution of the electron Y_e (thick line) and total lepton Y_l (thin) fractions as a function of central density for various models including GR. Left and right panels are without and with the NES effect, respectively. In the left panel, we see that the evolution of both Y_e and Y_l shows relatively good agreement among all models. Differences in both Y_e and Y_l are ∼0.01 at most. In addition, all numerical codes reproduce the neutrino trapping at nearly the same density of ρ ∼ 10¹² g cm⁻³ and Y_l stays nearly constant afterward.

Fig. 6. Central evolution of the electron (thick line) and lepton (thin line) fraction as a function of central density for various models. The left and right panels show this evolution without and with the NES effect being taken into account, respectively.

In the panel on the right, we see the role of NES. Because of the downscattering of neutrinos off the electrons, lower-energy neutrinos escape more efficiently and the lepton fraction Y_l becomes smaller compared to the model without NES. The parametrized deleptonization (PD) scheme enables Y_e to be nearly an identical value with the one in B95-SG. In G95-MSG, the upper envelope of Y_e shows a consistent value with the other models but with an oscillating behavior. Such oscillations are typical for simulations based on 12 energy groups. In corresponding spherically symmetric simulations (Mezzacappa & Bruenn 1993c; Liebendörfer et al. 2004) these oscillations are explained by an insufficient resolution of the Fermi-energy of the degenerate electron neutrinos in the late collapse phase, which does not affect the overall evolution of the dynamics. In the present 3D study, however, we did not have the resources to prove this using a corresponding simulation with 20 energy groups. In models G95-MSG and B95-SG, neutrino trapping is appropriately reproduced which can be seen from the constant Y_l evolution for ρ ≳ 10¹² g cm⁻³. The value of Y_l varies in models with PD scheme after neutrino trapping. This is because Y_ν is updated through the IDSA solver, but Y_e is taken from the PD formulation. Thus, the lepton number is not fully conserved with the PD scheme after neutrino trapping, but we consider that the overall picture is consistent with B95-SG/G95-MSG.

From the left panels of Fig. 1, the enclosed baryon mass with rest mass density larger than the neutrino trapping density (ρ ∼ 10¹² g cm⁻³) is larger (∼0.95 M_⊙) for models without NES (top-left panel) than for those (∼0.9 M_⊙) with NES (bottom-left panel). This can also be seen in the second column of panels, where we can see that the enclosed baryon mass of the central core is ∼0.05 M_⊙ larger for models without NES (∼0.65 M_⊙, top-right panel) compared to models with NES (∼0.60 M_⊙, bottom-right panel). This is because the additional pressure support from leptons leads to a more outward shock formation at bounce. It should, however, be noted that the impact of this scattering process depends also on the electron capture rate. Lentz et al. (2012b) reported that the down scattering can be blocked due to the production of lower-energy neutrinos when the up-to-date electron capture rate is used. Consequently, those models lead to the same results for the collapse phase, regardless of whether NES is included or not.

After bounce, we do not see any remarkable effects of NES. Although the shock evolution and heavy lepton-type neutrino luminosity show minor differences during the first ∼10 ms after bounce (see bottom two rows in Fig. 3) due to different shock position at bounce, the overall evolution shows no significant dependency on NES in time. We, however, note that the lack of visible effects of NES in the post-bounce phase might be due to the shortness of our simulation times. Just et al. (2018) recently reported, through their various models, that NES can turn a non-exploding model into a successful explosion. Nevertheless, such a noticeable effect appears at times (t_pb ≳ 300 ms) that are far beyond the end point of our simulations.

Calibration of parametrized deleptonization. Before more modern rates for electron captures can be implemented in the IDSA, the still missing implementation of NES in the IDSA would be important for the collapse phase. In combination with traditional electron capture rates, the inelastic scattering of neutrinos with electrons lowers the neutrino energy and accelerates deleptonization. In Lentz et al. (2012a), Pan et al. (2016), it has been shown that the central Y_e at core bounce is overestimated by about 10 − 15% if NES is not taken into account during the collapse. To effectively include NES, we follow the PD scheme described in Liebendörfer (2005), where the Y_e is parametrized by the baryon density and the chemical potential. Table 3 summarizes the parameters that we used in this paper. In Fig. 1, the PD scheme in X95-PG runs could lower the Y_e from 0.31 in X95-IG to 0.28, closely matching model G95-MSG, which includes NES.

In addition to electron fraction, neutrino stress is another important quantity that could affect the core size. The contribution of neutrino pressure is not only included in the IDSA solver, but also implemented in the PD scheme Liebendörfer 2005. The left panel of Fig. 7 shows a comparison of our 3D IDSA simulations with AGILE-BOLTZTRAN simulations without NES (orange lines). The nearly identical electron fraction profiles at core bounce indicate that the IDSA could capture the deleptonization correctly when the NES is off. Although the neutrino pressure contributes only a small fraction of the baryonic pressure, it could enlarge the inner core size by ∼30% (Fig. 7). Moreover, the right panel of Fig. 7 demonstrates that the PD scheme could mimic the NES effects and gives reasonable results compared to the AGILE-BOLTZTRAN simulation (the thick orange line; B95-SG). Furthermore, dashed lines present simulations of B95, F95-I, and F95-PG without the contribution from neutrino stress (with label “N”), leading to an inner core that is smaller by ∼30%, consistent with the expected value.

Fig. 7. Influence of the neutrino stress in the Ye profiles at bounce. When neutrino stress is not included (dashed lines, i.e., N-labeled models) the inner core is ∼30% smaller. Left panel: IDSA Newtonian calculations, right panel: GR calculations including neutrino-electron scattering.

4.3. The role of general relativity

General relativity plays an essential role in CCSNe. The most detailed 1D simulations to date, in which the GR Boltzmann transport equation was solved in spherical symmetry, were performed by Wilson (1971), Mezzacappa & Matzner (1989), Yamada (1997), Liebendörfer et al. (2001), Sumiyoshi et al. (2005), Lentz et al. (2012a), and their 3D counterparts, with the multi-angle and multi-energy neutrino transport, by Sumiyoshi & Yamada (2012). Although it is not a full Boltzmann transport scheme, the variable Eddington factor (VEF) method has been coupled to a full GR hydro solver (e.g., Müller et al. 2012, and references therein). In this method, one can self-consistently determine the closure relation from a model Boltzmann equation that is integrated to iteratively obtain the solution up to the higher moments (i.e., the Eddington tensor) until the system converges.

Bruenn et al. (2001), Liebendörfer et al. (2001), with baseline neutrino opacities, and Lentz et al. (2012a), with the currently best available ones, presented evidence that the average neutrino energy of any neutrino flavor during the shock reheating phase increases when switching from Newtonian to GR hydrodynamics. They also pointed out that the increase is larger in magnitude compared to the decrease due to redshift effects and gravitational time dilation. In addition, Lentz et al. (2012a) showed that the omission of observer corrections in the transport equation lessens the triggering of neutrino-driven explosions. In these fully fledged 1D simulations, a commonly observed disadvantageous aspect of using GR to trigger neutrino-driven explosions is that the residency time of material in the gain region becomes shorter due to the stronger gravitational pull. As a result of these competing ingredients, in the end GR works against the triggering of the neutrino-driven explosions in 1D. In fact, the maximum shock extent in the post-bounce phase is shown to be ∼20% smaller when switching from Newtonian to GR hydrodynamics (Bruenn et al. 2001; Liebendörfer et al. 2001; Sumiyoshi et al. 2005; Lentz et al. 2012a).

From Fig. 3, we can see that the mean PNS radii are ∼10% smaller in (effective-)GR models compared to Newtonian ones (see, e.g., the second and fourth panels from the top in the leftmost column), confirming the same effect also in 3D simulations. GR also affects the shock evolution. Irrespective of NES, the shock front appears ∼10 km more inward in (effective-)GR models. Furthermore, because of the more compact and hotter PNS, neutrino luminosities in all flavors become ∼20 % higher. AGILE-BOLTZTRAN takes into account both the gravitational redshift and Doppler shift terms. Therefore, we can explore the influence of GR on the emergent neutrino energies by comparing B95-S and B95-SG. As a consequence of the competition between the gravitational red-, Doppler-shift, and hotter PNS surface, the mean neutrino energies are higher in B95-SG than in B95-S by ∼0.5 MeV. This simply reflects the effect of GR, leading to a more compact PNS core, which can also be confirmed in Fig. 4. This figure shows a comparison of the central density evolution between Newtonian (left) and GR (right) models. In GR models, the central density shows values that are systematically higher by ∼20 − 30%. To see if these GR effects work advantageously on the shock revival, we need to perform simulations that cover more of the post-bounce phase.

4.4. The role of rotation

With the advent of multidimensional SN simulations the role of rotation is gaining increased attention from the community (Marek & Janka 2009; Suwa et al. 2010; Yokozawa et al. 2015; Takiwaki et al. 2016; Janka et al. 2016; Gilkis 2018; Summa et al. 2018), as it can have sizable effects on the emergence of an explosion as well as on important SN observables. For example, gravitational waves produced from the development of triaxial instabilities triggered by rotation will carry information about the supernuclear EOS that is encoded in their waveform. From these, we could impose constraints on rotation rates of the iron core and the PNS, which are not well known.

Following the setup presented in Yokozawa et al. (2015) we implemented a shellular rotation (Eriguchi & Mueller 1985) assigning an angular velocity Ω at each cylindrical radius r in the form of

(31)

where Ω₀ is the angular velocity at the center and r₀ is a characteristic distance that controls the amount of differential rotation. For the limit of small values of r₀, Eq. (31) approaches a configuration with constant specific angular momentum, while for large values it tends to a rigidly rotating model. We chose r₀ = 1000 km and Ω₀ = π rad s⁻¹. The former ensures a nearly rigid rotation of the core of the star, while the latter corresponds to a spin period of 2 s. We note that this is roughly a factor 50 more rapid than the predictions by detailed stellar evolution models (Heger et al. 2005) that include angular momentum transfer by magnetic processes. Nevertheless, many publications employ such a rapid rotation profile, even if it is probably only compatible with a differentially rotating core in hypernovae (Burrows et al. 2007) and/or a fast rotating black hole at the origin of gamma-ray bursts (MacFadyen & Woosley 1999). In order to be able to compare with the literature we kept the same rotational profile. We also highlight here that Eq. (31) generates a rotational profile in cylindrical shells, which was not applicable in ELEPHANT. As explained in Sect. 2.1, this code simulates the inner region of the star with a 3D Cartesian mesh while its boundary conditions are provided by a simultaneous calculation with the 1D code AGILE-IDSA. This leads us to impose only spherical rotational profiles for the 3D inner region. Therefore, ELEPHANT applies Eq. (31) with r and r₀ representing spherical radii instead of cylindrical radii. As a consequence, the initial rotational profile of ELEPHANT is different from those of SPHYNX and FLASH as it rotates more slowly (∼6%).

Figure 8 presents an overview of the average shock position and PNS radius (left column), electron-neutrino luminosity (second column), heavy-lepton neutrino luminosity (third column), and mean energies for all neutrino species (last column) for the calculated rotating models without (IR) and with (PR) parametrized deleptonization. In order to clearly see the influence of rotation, we added a third row to Fig. 8 with the results of X07-P models, which include the same physics as the rotating models but without rotation. Overall, there is good agreement among the codes for the rotating models. All three do not show signs of a rapidly stalled shock within the simulated time, in contrast to the non-rotating models where the shock halts at about 100–150 km within the first 20 ms. As expected, rotation helps the expansion of the shock due to the deposition of angular momentum in outer layers (Fryer & Warren 2004; Suwa et al. 2010; Nakamura et al. 2014). There is good agreement in the neutrino and anti-neutrino luminosities among all thee codes, although SPHYNX tends to be slightly less luminous than the other two codes for anti-neutrinos (second column), which correlates with the lower mean energy deposited in that neutrino species (fourth column). Nevertheless, the main difference among codes can be seen in the heavy-lepton neutrino luminosity (third column), where SPHYNX and FLASH show good agreement, but ELEPHANT produces ∼40% more luminosity. The reason for this is the different rotational profile of ELEPHANT, which produces a more compact PNS, consistent with its central density evolution shown in Fig. 9. As expected, ELEPHANT simulations tend to form more compact objects than those of FLASH. SPHYNX obtains intermediate results, where the biggest differences are ∼15%.

Fig. 8. Similar to Fig. 3 but for rotating models X07-IR (first row) and X07-PR (second row). Non-rotating models X07-P (third row) are shown here for comparison purposes.

Fig. 9. Central density evolution for rotating models without (left panel) and with (right panel) parametrized deleptonization.

The angular momentum transfer is crucial in rotating scenarios. We note that the angular momentum transport is strongly enhanced after bounce, where the shock launch aids the transfer of fast-rotating matter to larger radii, driving a rapid decrease of the angular momentum of the inner core, as it is transported to the surface of the neutron star. This mechanism excites the m = 1 spiral-arm modes, leaving a signal in the fluid that we can see in Figs. 10 (for IR models). In this latter figure, at t_pb = 50 ms, we plot the deviation from the spherical average of the density (left), pressure (center), and radial velocity (right) in a thin slice of the 3D domain along the equatorial plane. We can see in all three codes that a spiral, high-density, high-pressure wave is moving outwards from the center of the system. SPHYNX and FLASH produce very similar results, both in magnitude and size of the spiral perturbation at the same time. Due to the different rotational profile, ELEPHANT results are not directly comparable to those of the other two codes. Nevertheless, a spiral shape is still visible on all three magnitudes, albeit ∼15% more extended. Our results are similar to those of Takiwaki et al. (2016; see their Fig. 3).

Fig. 10. Deviations from the spherical average of density (left panels), pressure (center panels), and radial velocity (right panels) profiles for rotational models X07-IR at 50 ms.

Additionally, we can also follow the transfer of angular momentum directly during the whole evolution, plotting the difference of specific angular momentum with respect the initial specific angular momentum distribution as a function of the enclosed mass. This is represented in Fig. 11 for ELEPHANT, FLASH, and SPHYNX (from top to bottom) at different stages of the collapse (left column) and at different times after bounce (right column). If we focus on the first column, we can see that even during collapse, some angular momentum is lost in the inner core for all three codes. This may be due to two factors that occur at the same time: the appearance of shear and the lack of angular momentum conservation by the code. Some shear might be expected during the collapse of a rotating star due to the stronger centrifugal forces experienced by the material with higher angular momentum (i.e., in the equatorial plane). As a consequence, its collapse proceeds slower than material with lower angular momentum (e.g., in the polar regions) and a certain quantity of shear appears, causing angular momentum to pile up at larger radii (this is especially visible at t_pb = 0 ms in the SPHYNX run, in the small bump from M ∼ 0.8 M_⊙ up to ∼1.3 M_⊙ in Fig. 11). Additionally, some angular momentum is lost due to numerical diffusion. This is more relevant in the mesh-based codes ELEPHANT and FLASH, while SPHYNX does conserve it better. This is an expected behavior of an SPH code, which, being purely Lagrangian, is specially constructed to conserve momentum and energy. It is worth noting though, that some angular momentum is also lost in SPHYNX due to the calculation of the gravitational force, which, to some extent, always breaks the perfect conservation properties of the SPH formalism. We can estimate the total angular momentum that is lost due to numerics and shear by integrating the losses and gains of the specific angular momentum as separate functions of the enclosed mass. If angular momentum is perfectly conserved, both magnitudes should be equal and therefore their difference should be zero. We found that, up to bounce, SPHYNX has a losses-to-gain difference below 1%, while the value for FLASH is below 3%, and the value for ELEPHANT is below 4%. This percentage gives an upper limit to the angular momentum losses at the (low) spatial resolution at which these simulations have been performed. Increasing the resolution, this ratio improves considerably for all three codes. Unfortunately we cannot clearly disentangle the undesired angular momentum loss on the scale of the dimensions of the PNS from unavoidable angular momentum loss due to dissipation of turbulence on smaller scales. But in our investigation of the collapse phase and bounce we consider the loss-to-gain difference as a reasonable estimate for undesired angular momentum losses because our simulations do not develop turbulence during the very short collapse of our perfectly spherically symmetric progenitor models.

Fig. 11. Specific angular momentum differences with respect to the initial specific angular momentum distribution at different central densities during collapse (left panels) and at different post-bounce times (right panels). From top to bottom panels, the represented simulations correspond to ELEPHANT, FLASH, and SPHYNX codes.

After bounce (second column) it is clear that angular momentum is extracted from the inner core and transported to the outer layers. The time and the regions at which angular momentum changes differ among codes. The presence of convection is the reason for the complex structure of these profiles. The different treatments of dissipative terms in the different codes and the low spatial resolution may also have an important impact. Despite these quantitative differences, it can be seen in all three codes how angular momentum moves up along the enclosed mass axis, and, in fact, reaches similar enclosed masses, that is, ∼1.4 M_⊙, in FLASH and SPHYNX at about the same time (t ∼ 50 ms).

4.5. Computational resources

Using different codes with different levels of approximation implies great differences in computational cost, which, in turn, affect the spatial resolution that we can afford. In Table 4 we provide some details regarding the CPUh consumption, parallelization methods, number of fluid elements, and resolution of each code in this work in order to give some idea of the computational costs. Some of the calculations were done in different facilities, so the used resources presented in Table 4 are merely indicative. Additionally, the computational efficiency of each code is obviously highly dependent on its parallelization methods and its implementation. Hence, the values presented in Table 4 reflect the current state of development of these codes.

ELEPHANT provides (highly efficient) simulations with constant spatial resolution (1 km), which is very good at the shock position but coarse at the PNS. This is mitigated in FLASH via the usage of adaptive mesh refinement (AMR). The resolution at the PNS is greatly improved and a third-order scheme is used at the price of increased computational resources. fGR1 uses nested meshes of finer resolution, so there is no AMR overhead. Nevertheless, the inclusion of full GR and M1 treatments requires considerable computational resources. On the other hand, SPHYNX reaches similar resolution in the PNS with considerably less fluid elements, but, as in SPH codes resolution follows density, more diluted regions like that of the shock are much coarser. In the case of ELEPHANT and SPHYNX, the lack of resolution in some regions of the core collapse scenario can be compensated with higher amounts of fluid elements. Subsequently, the viability of the simulation relies on its computational efficiency and the scaling of the code in a parallel implementation. It is clear that the use of one code over another can be beneficial depending on the focus of the study. For example, increasing the amount of SPH particles up to more “production-simulation-like” amounts (e.g., 2 × 10⁷) will increase the spatial resolution by a factor ∼4.6 reaching ∼100 m resolution in the PNS. On the other hand, ELEPHANT can provide a very high resolution (∼500 m) in the gain region and at the same time include magnetic fields.

5. Conclusions

In this work, we studied the same core collapse scenario in 3D with four different hydrodynamics codes including, for the first time, Eulerian and Lagrangian codes in the same study. In addition, we also compared the four hydrodynamics codes with a spherically symmetric code with Boltzmann transport. Within this framework, we varied three parameters in order to explore their different impacts on different hydrodynamics codes, namely the neutrino-electron scattering, a monopole GR correction, and rotation. Additionally, we discussed the differences between the approximate treatments and more detailed calculations, that is, IDSA+Parametrized deleptonization versus M1 neutrino transport, and the monopolar GR correction versus full GR.

The comparison of different neutrino transport methods in a SN simulation revealed that all our codes are able to reproduce a reasonable overall evolution of the early post-bounce phase in 3D. The strong coupling and feedback between the PNS compactness, shock dynamics, neutrino luminosities, and mean energies is in all cases consistent. Despite their different discretization methods of the Euler equations, including radically different domain decompositions, parallelization approaches, and the use of static or nested meshes, AMR or SPH, the deviations of the results are of the order of 10%, and sometimes 20%. This comparison leads us to propose a requirement for further investigation of the luminosities of antineutrinos in the IDSA, which tend to be rather high. Furthermore, the leakage scheme for the heavy neutrinos may not accurately represent the corresponding energy loss under all conditions of the PNS. Here, a better scheme would be desirable (e.g., Perego et al. 2014; Takiwaki et al. 2014). The most computationally expensive full GR code still seems to struggle with resolution and performance issues. The whole ensemble of codes however presents many virtues if referred to in combination: excellent resolution of the PNS and angular momentum conservation (SPHYNX), flexibility and open accessibility for efficient 3D SN simulations with GPU-acceleration (FLASH), superior shock resolution and efficient inclusion of magnetic fields (ELEPHANT), fully GR hydrodynamics and neutrino transport (fGR1). The results presented here verify the usefulness of these four codes with respect to their coupling with the neutrino treatment and to the implementation of gravity as previously discussed.

We also checked the influence of neutrino-electron scattering, which plays a significant role during collapse phase. Whether we implement the inelastic scattering kernel explicitly in the source terms or control the deleptonization process in a parametrized way, we confirm that the lepton profiles within the inner core show converged results among all models. Furthermore, the hydrodynamics and neutrino profiles post-bounce show consistent evolution with deviations of ∼10%, which is an acceptable level considering that we use totally independent codes.

The effects of GR are investigated by comparing the Newtonian and (effective-)GR models. Due to the stronger gravitational force, we see more compact PNS cores and shock positions in GR models. This is consistent with previous studies (Bruenn et al. 2001; Liebendörfer et al. 2001; Sumiyoshi et al. 2005; Lentz et al. 2012a). From hotter and more compact PNS cores, the emergent neutrinos show systematically higher neutrino energies and luminosities in all flavors, even after taking into account the effect of gravitational redshift (e.g., comparing B95-S and B95-SG). Within our simulation times, we do not see any remarkable difference between full- and effective-GR models, and therefore we consider the usage of an effective-GR potential as justified. However, our simulation time is short and longer simulations are required to see for how far the effective-GR potential works properly and to determine the time when those two methods diverge. Furthermore, it would also be interesting to check under which conditions GR acts advantageously on the shock revival compared to the Newtonian case.

The inclusion of rotation helps the expansion of the shock in general. Even when using different rotation profiles (ELEPHANT vs. FLASH and SPHYNX), the shock position at 50 ms can be 26 − 70% further than in models without rotation. We confirmed the presence of a m = 1 spiral perturbation that drives matter from the inner core into the gain region, aiding the explosion. This is partially compensated by a less dense and cooler PNS, but the net balance between both effects favors the expansion of the sock. Additionally, we could evaluate the performance of the different hydrodynamics codes with respect to the angular momentum transport. With the present (low) resolution in this study the angular momentum conservation was better than 4% .

Additionally, we were able to clearly demonstrate that an SPH code with a spectral neutrino treatment produces results comparable to those of Eulerian codes. Interestingly, this opens the possibility to consider SPH as a competitive, efficient method to simulate core-collapse SNe with adaptive resolution and open boundaries. This is particularly useful when one aims to study the PNS at very high resolution, when angular momentum conservation is capital, or when the PNS kick is the focus of an investigation. We note, that there is remarkable agreement between the results of SPHYNX and FLASH in many simulated models. In particular when the effective GR potential is used. Both codes provide a high spatial resolution at the PNS, which is important for an accurate compactness of the PNS and thus the overall evolution of the system.

It should be noticed that there are several multi-dimensional effects that deserve further investigation, such as turbulence, standing accretion shock instability, magnetic fields, and gravitational wave emission.

Additionally, extending this kind of comparison to other progenitors and longer timescales is of utmost importance. However, this kind of simulation is extremely costly in terms of computational resources and as a consequence has been relegated. We think, though, that the time to perform detailed comparisons among complex production codes and across many different simulations is getting closer. It is in this spirit that our work provides the first outlook on SN code comparisons in 3D.

Acknowledgments

The authors want to thank S. C. Whitehouse for his extensive work on ELEPHANT. This work has been supported by the European Research Council (FP7) under ERC Advanced Grant Agreement No. 321263 - FISH, by the Swiss Platform for Advanced Scientific Computing (PASC) project DIAPHANE (RC and KCP) and SPH-EXA (RC), by ERC Starting Grant EUROPIUM-677912 (TK), and by JSPS KAKENHI Grant Number JP17H01130 and JP17H06364 (TK). A part of numerical computations were carried out on Cray XC30 at CfCA, National Astronomical Observatory of Japan (TK). The authors acknowledge the support of sciCORE (http://scicore.unibas.ch/) scientific computing core facility at University of Basel, the Swiss National Supercomputing Center CSCS (http://cscs.ch) under grant No. 661 and No. 840, where these calculations were performed. FLASH was in part developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago. Analysis and visualization of the FLASH simulation data were completed using the yt analysis toolkit (Turk et al. 2011).

References

All Tables

All Figures

	Fig. 1. Rest mass density (left column) and electron/lepton fraction profiles as a function of the enclosed baryon mass for various models at bounce. Upper and lower rows are without and with taking into account the NES effect, respectively.
In the text

	Fig. 2. Time evolution of mass accretion rates measured at R = 200 km. Different colors represent simulations with different codes. Each column shows a set of simulations (see Table 2).
In the text

	Fig. 3. Time evolution of averaged shock radius, PNS radius, neutrino luminosities, and mean energies for different runs. Each row describes a comparison of different physics input as described in Table 2. Different color represents a simulation with a different hydrodynamics code.
In the text

	Fig. 4. Central density evolution in all codes within a Newtonian (left panel) and GR (right panel) framework. The dashed line shows AGILE-BOLTZTRAN data before the extrapolation to the center (cf. text).
In the text

	Fig. 5. Central entropy as a function of central density (left panel) and as a function of time (right panel) for GR models including NES.
In the text

	Fig. 6. Central evolution of the electron (thick line) and lepton (thin line) fraction as a function of central density for various models. The left and right panels show this evolution without and with the NES effect being taken into account, respectively.
In the text

	Fig. 7. Influence of the neutrino stress in the Ye profiles at bounce. When neutrino stress is not included (dashed lines, i.e., N-labeled models) the inner core is ∼30% smaller. Left panel: IDSA Newtonian calculations, right panel: GR calculations including neutrino-electron scattering.
In the text

	Fig. 8. Similar to Fig. 3 but for rotating models X07-IR (first row) and X07-PR (second row). Non-rotating models X07-P (third row) are shown here for comparison purposes.
In the text

	Fig. 9. Central density evolution for rotating models without (left panel) and with (right panel) parametrized deleptonization.
In the text

	Fig. 10. Deviations from the spherical average of density (left panels), pressure (center panels), and radial velocity (right panels) profiles for rotational models X07-IR at 50 ms.
In the text

	Fig. 11. Specific angular momentum differences with respect to the initial specific angular momentum distribution at different central densities during collapse (left panels) and at different post-bounce times (right panels). From top to bottom panels, the represented simulations correspond to ELEPHANT, FLASH, and SPHYNX codes.
In the text