Core-collapse supernova simulations in one and two dimensions: comparison of codes and approximations
We present spherically symmetric (1D) and axisymmetric (2D) supernova simulations for a convection-dominated 9 and a 20 progenitor that develops violent activity by the standing-accretion-shock instability (SASI). We compare in detail the Aenus-Alcar code, which uses fully multidimensional two-moment neutrino transport with an M1 closure, with a ray-by-ray-plus (RbR+) version of this code and with the Prometheus-Vertex code that employs RbR+ two-moment transport with a Boltzmann closure. Besides testing consequences of ignored non-radial neutrino-flux components in the RbR+ approximation, we also discuss the influence of various transport ingredients applied or not applied in recent literature, namely simplified neutrino-pair processes, neutrino-electron scattering, velocity-dependent and gravitational-redshift terms, and strangeness and many-body corrections for neutrino-nucleon scattering. Alcar and Vertex show excellent agreement in 1D and 2D despite a slightly but systematically smaller radius (1 km) and stronger convection of the proto-neutron star with Alcar. As found previously, the RbR+ approximation is conducive to explosions, but much less severely in the convection-dominated 9 case than in the marginally exploding 20 model, where the onset time of explosion also exhibits big stochastic variations, and the RbR+ approximation has no distinctly stronger supportive effect than simplified pair processes or strangeness and many-body corrections. Neglecting neutrino-electron scattering has clearly unfavorable effects for explosions, while ignoring velocity and gravitational-redshift effects can both promote or delay the explosion. The ratio of advection timescale to neutrino-heating timescale in 1D simulations is a sensitive indicator of the influence of physics ingredients on explosions also in multidimensional simulations.
keywords:hydrodynamics – instabilities – radiative transfer – supernovae: general – neutrinos
Rarely any astrophysical phenomenon has been studied as long and as intensely using the most powerful available supercomputers as the explosion mechanism of core-collapse supernovae (CCSNe; see, e.g., Foglizzo et al., 2015; Janka et al., 2016; Müller, 2016, for recent reviews). After several decades of research the neutrino-driven mechanism, first suggested by Colgate & White (1966) and later in its modern version by Bethe & Wilson (1985), remains the leading candidate responsible for the revival of the initially stalled shock wave in ordinary CCSNe. Although current three-dimensional simulations lend support to the scenario that neutrino heating powers the explosion (Takiwaki et al., 2014; Melson et al., 2015a; Melson et al., 2015b; Lentz et al., 2015; Roberts et al., 2016; Ott et al., 2017) important questions remain to be studied in more detail and are vividly investigated by a large number of researchers worldwide. These questions concern, e.g., the relative importance of the standing accretion-shock instability (SASI; Blondin et al., 2003) and post-shock convection (e.g. Fernández, 2015), the impact of improved neutrino interactions (e.g. Kotake et al., 2018), the role of turbulence (e.g. Radice et al., 2018), the dependence on the progenitor (e.g. Sukhbold et al., 2016) and its pre-collapse asymmetries (e.g. Müller et al., 2017), the quest for a universal explosion criterion (e.g. Ertl et al., 2016), the information contained in the emitted neutrinos (e.g. Tamborra et al., 2017), gravitational waves (e.g. Richers et al., 2017a), and remnant structure (e.g. Ono et al., 2013), or the role of rotation and magnetic fields (e.g. Obergaulinger & Aloy, 2017).
All these research lines depend in some way or another on results obtained from time-dependent models of the supernova core, which are provided by multidimensional (neutrino-) hydrodynamics simulations. One of the major obstacles for detailed CCSN models is the neutrino transport. Ideally, the latter requires to solve along with the three-dimensional hydrodynamics equations a Boltzmann equation depending on six phase-space coordinates. However, since this is not possible with present supercomputer capabilities, at least not for reasonable resolution and number of time steps, various schemes and approximations have been proposed in the past, which differ considerably concerning their computational efficiency and accuracy.
The computationally most efficient methods of handling neutrino processes are light-bulb (e.g. Hanke et al., 2012; Couch & Ott, 2015) and leakage (Ruffert et al., 1996; O’Connor & Ott, 2010; Perego et al., 2016) schemes, which add almost no cost to the hydrodynamics equations since the neutrino source terms entering the latter are estimated based on the hydro-/thermodynamic information only. More involved schemes, such as the fast multigroup transport (Müller & Janka, 2015) or the one designed by Scheck et al. (2008), solve the conservation equation for neutrino energy assuming a pre-defined profile for the ratio of the flux density to energy density. Another scheme with relatively high computational efficiency is the isotropic diffusion source approximation (IDSA; Liebendörfer et al., 2009; Takiwaki et al., 2014; Suwa et al., 2016) that decomposes the neutrino distribution into a trapped and a free-streaming component, which are separately evolved using a parabolic diffusion equation and an elliptic equation, respectively. One of the most often used schemes, also outside of CCSN theory, is flux-limited diffusion (FLD; Levermore & Pomraning, 1981; Dessart et al., 2006; Bruenn et al., 2016; Dolence et al., 2015), which evolves the energy density, , of neutrinos assuming the flux density, , to be given locally by a generalized (“flux-limited”) diffusion law that limits the flux to remain causal, ( being the speed of light).
More recently, the so-called M1 scheme (Minerbo, 1978; Levermore, 1984; Pons et al., 2000; Audit et al., 2002; Shibata et al., 2011) was implemented in a number of neutrino-transport codes (O’Connor, 2015; Foucart et al., 2015; Sekiguchi et al., 2015; Just et al., 2015; Kuroda et al., 2016; Skinner et al., 2016). Being closely related to FLD, M1 evolves additionally to the latter the flux densities. The resulting two-moment system – energy and flux density are the 0th- and 1st-order angular moments of the specific neutrino intensity – is augmented with a local closure relation expressing the 2nd-order moments (i.e. the neutrino pressure tensor) as a function of the evolved moments. The M1 scheme comes with the convenient property of being hyperbolic, which, together with the fact that fluid velocities are high in CCSNe, allows to advance the equations using explicit time stepping. This, in turn, makes it computationally feasible to solve the unconstrained, fully multidimensional transport equations.
Finally, in contrast to the aforementioned approximate schemes, more involved but also considerably more expensive schemes, so-called Boltzmann solvers, attempt to resolve the full angular dependence of the radiation field. This can be achieved, e.g. by employing tangent ray (Burrows et al., 2000; Rampp & Janka, 2002; Buras et al., 2006), discrete ordinate (Mezzacappa & Bruenn, 1993; Yamada et al., 1999; Liebendörfer et al., 2004; Nagakura et al., 2018), or Monte Carlo methods (Janka & Hillebrandt, 1989; Richers et al., 2017b; Foucart, 2018) to evolve the Boltzmann equation either directly or coupled to a two-moment system.
One measure to make expensive solvers more affordable in axisymmetric (2D; e.g. Buras et al., 2006; Müller et al., 2012a; Bruenn et al., 2016; Nakamura et al., 2015) and three-dimensional (3D; e.g. Takiwaki et al., 2014; Melson et al., 2015a; Lentz et al., 2015; Müller et al., 2017; Wongwathanarat et al., 2017) simulations is to employ the ray-by-ray(-plus) approximation (RbR+), which neglects the evolution of the non-radial flux components and therefore allows to evolve basically one-dimensional (1D) transport problems on each angular bin along the radial directions in a quasi-decoupled manner. This leads to a major boost of computational performance, particularly for a scheme where the time integration requires the inversion of matrices spanned over the entire grid, because then the total number of operations is reduced by much more than a factor of three and the code can be efficiently parallelized. The ray-by-ray approach is naturally motivated by the rather spherical geometry of the SN core and resulting subdominance of non-radial neutrino fluxes, and by the notion that local anisotropies enhanced by RbR+ could average out and thus only leave a small imprint on the angle- and time-averaged dynamics (Buras et al., 2006). On the other hand, the suspicion was raised (Dolence et al., 2015; Sumiyoshi et al., 2015) that the RbR+ approximation could enhance the feedback between the fluid flow and neutrino radiation, which could induce explosions more readily. However, so far the RbR+ approximation was investigated self-consistently, i.e. by comparing time-dependent simulations with and without using RbR+, only using one code (Skinner et al., 2016) and more detailed exploration is warranted.
In addition to the large number of available numerical approximations of the neutrino transport, probably an even larger number of (micro-)physics ingredients exists in various degrees of sophistication. For instance, while some studies use a large set of state-of-the-art neutrino interactions (e.g. Melson et al., 2015a; Lentz et al., 2015; Kotake et al., 2018), many other studies rely on a smaller and more basic set of interactions (often adopted from Bruenn, 1985) and neglect or simplify numerically cumbersome reactions that couple different neutrino energies, such as inelastic scattering of neutrinos off electrons and positrons as well as pair processes (e.g. electron-positron annihilation and nucleon-nucleon bremsstrahlung). Similarly, frame-dependent effects such as Doppler- and gravitational energy-shifts are sometimes neglected for numerical convenience by dropping all velocity-dependent terms and energy derivatives in the evolution equations. Other features, such as strangeness (Horowitz, 2002) or many-body (Horowitz et al., 2017) corrections to the neutral-current cross sections for neutrino-nucleon scattering can be implemented rather easily and have recently found their way into a number of simulations (e.g. Melson et al., 2015a; Bollig et al., 2017; Burrows et al., 2018; Kotake et al., 2018).
With the growing number of numerical models, also the demand grows to perform systematic comparisons between codes, methods and approximations. Method comparisons, such as the ones conducted by Liebendörfer et al. (2005); Lentz et al. (2012b, a); Müller et al. (2012a); Richers et al. (2017b); O’Connor et al. (2018); Pan et al. (2018); Cabezón et al. (2018), might not only help improving the employed algorithms on all sides, e.g. by unmasking previously unknown deficiencies. They may also prove useful to other groups in locating the reasons why their codes give different results. Moreover, they may provide instrumental help in developing and gauging code extensions or new codes. Unfortunately, going beyond spherically symmetric models when comparing time-dependent CCSN simulations is significantly more challenging because of the appearance of (numerical and physical) perturbations, fluid instabilities, turbulence, and stochasticity. The last mentioned issue, the importance of which was emphasized by Cardall & Budiardja (2015), may be a particular cause of complication if simulations are heavily time consuming and expensive and therefore only allow for a small number of runs to be conducted.
In this paper we use the M1-based code Aenus-Alcar (Just et al., 2015; Obergaulinger, 2008) to compare in spherical- and axisymmetry against the well-known code Prometheus-Vertex (Rampp & Janka, 2002; Buras et al., 2006) that employs the RbR+ approximation with a Boltzmann closure, and to test the ramifications of the RbR+ approximation in axisymmetry together with those of other frequently used modeling simplifications in the transport sector. More specifically, we address the following questions:
How good is the agreement between the M1 code Alcar and the Boltzmann solver Vertex in 1D given exactly the same input physics?
How good is the corresponding agreement in 2D?
What is the impact of RbR+ on the explodability and on proto-neutron star (PNS) convection, how does this impact depend on the stellar progenitor, and how does it compare to the other considered modeling variations?
What difference does it make if neutrino-electron scattering is neglected, pair processes are simplified, and strangeness as well as many-body corrections are included?
What happens if velocity-dependent terms as well as gravitational redshift are neglected in the transport?
What is the impact of stochasticity, e.g. how large is the scatter of explosion times if simulations are repeated several times?
2 Numerical methods and model setup
2.1 Governing equations and discretization schemes
In the following we outline the main features of the employed simulation tools, Alcar and Vertex. For in-depth descriptions of these codes we refer the reader to Just et al. (2015) as well as Rampp & Janka (2002) and Buras et al. (2006), respectively.
Both codes employ a Godunov-type finite-volume scheme in spherical polar coordinates to solve the equations of Newtonian hydrodynamics with an effective general relativistic gravitational potential (Marek et al., 2006). The equations read (with vector indices running over radial, , polar, , and azimuthal, , coordinate, and ):
where , and are the baryon density, fluid velocity, electron fraction, total (kinetic plus internal) gas energy density, gas pressure, and gravitational potential, respectively, and the neutrino-related source terms are given by (with neutrino-energy coordinate and baryonic mass constant )
in terms of the lapse function and the 0th- and 1st-order angular moments of the collision integral measured in the comoving (i.e. fluid rest) frame, and , respectively, where the index runs over all six neutrino species.
Both and are computed from radial profiles of angle averaged fluid- and neutrino-related quantities using the prescriptions of Marek et al. (2006) (for the case “A” potential) and Rampp & Janka (2002), respectively. We do not include spherical harmonics terms higher than the monopole in the gravitational potential, i.e. gravity is considered to be spherically symmetric.
The hydrodynamics module coupled to Alcar, Aenus (Obergaulinger, 2008), offers several choices for spatial reconstruction methods, Riemann solvers, and time integration schemes. In this study, we employ spatial reconstruction as adapted from the piecewise-parabolic method (PPM; Colella & Woodward, 1984) in the version by Mignone (2014) that retains high-order accuracy near the coordinate center and polar axis111More specifically, we use the cell-centered version “PPM” of Mignone (2014) to reconstruct the primitive variables and employ flattening near strong shocks as in Colella & Woodward (1984) but no steepening near contact discontinuities. We do not perform characteristic tracing as in the original PPM method., the HLLC Riemann solver everywhere except near coordinate-aligned shocks where we switch to the HLLE solver (e.g. Toro, 1997), and a dimensionally unsplit 2nd-order Runge-Kutta time-integration scheme.
The neutrino transport in both codes is based on the multi-group evolution of the specific (i.e. per unit of neutrino energy) energy density, , and specific flux density, , both measured in the comoving frame of the fluid, for the three species electron-neutrinos , electron-antineutrinos , as well as representative of the four remaining heavy-lepton neutrinos. The latter is governed by (suppressing indices ):
where are the 2nd- and 3rd-order angular moment tensors, respectively. Compared to the purely Newtonian counterparts, Eqs. (3) contain corrections for general relativistic redshift and time dilation. We use the same corrections in Alcar as those in Vertex, which have been motivated in Rampp & Janka (2002)222Specifically, equations (3) can be recovered from the general relativistic moment equations (e.g. Shibata et al., 2011; Cardall et al., 2013, where necessary after transforming the moments into the comoving frame) by assuming the only general relativistic metric component to be the lapse function, , performing the replacement , and dropping all terms proportional to and . and shown, e.g. in Liebendörfer et al. (2005); Müller et al. (2012a), to lead to results that agree reasonably well with fully general relativistic results.
In Alcar, the higher-order moments and needed to close the set of moment equations are expressed as functions of the moments by
where , , and the quantities and are determined by the chosen one-dimensional closure. For the latter, we use the Minerbo-closure and refer to Minerbo (1978) and Just et al. (2015) for explicit expressions of and , respectively. The resulting hyperbolic system of equations is solved in close analogy to the hydrodynamics equations, namely on the same grid and with the same reconstruction scheme, and an HLL Riemann solver. The time integration of the potentially stiff interaction source terms, however, is done in a mixed explicit-implicit manner in order to ensure numerical stability (see Appendix A for details).
Vertex solves the same set of equations, Eqs. (1) and (3), for hydrodynamics and neutrino transport, respectively, except for the following differences: Vertex obtains the higher-order moments using a variable-Eddington-factor approach, i.e. by solving a (slightly simplified) Boltzmann equation in addition to the moment equations. This approach is comparable in accuracy to solving the Boltzmann equation directly for a spherically symmetric stellar background and is then constrained to be used in the RbR+ mode. The latter assumes that the radiation field is axisymmetric around each radial ray, leading to vanishing non-radial components of the flux density vector, , while lateral advection of neutrinos by the fluid as well as lateral neutrino-pressure forces are still taken into account (see Buras et al., 2006, for details). Vertex evolves Eqs. (3) using finite-difference methods and fully implicit time integration. For the hydrodynamics part the well-known Prometheus code (Fryxell et al., 1989) is applied, which employs the original PPM method (Colella & Woodward, 1984) to integrate the hydrodynamics equations with 2nd-order Strang-type dimensional splitting.
2.2 Model setup and neutrino interaction channels
We consider two stellar progenitor models in this study, one with rather high zero-age main sequence (ZAMS) mass of , model s20 (Woosley & Heger, 2007), and one with lower ZAMS mass of (model 9.0A of Woosley & Heger, 2015333We use the slightly upgraded version with respect to Woosley & Heger (2015) that was also employed by Sukhbold et al., 2016 and that is available at https://wwwmpa.mpa-garching.mpg.de/ccsnarchive/data/SEWBJ_2015/index.html.), which we denote here as model s9. The simulations are initialized using the density, temperature, and electron fraction from the progenitor data. The 2D models are set up by mapping from 1D simulations at around ms after bounce and adding random perturbations to the density, , in each grid cell. In doing so, we ignore recent findings (e.g. Couch et al., 2015; Müller et al., 2016) that rather strong perturbations with low angular order can be present before collapse, which may help initiating the shock runaway. Since in this study we are only concerned about the impact of modeling variations on the post-bounce evolution, all two-dimensional models are mapped from the same 1D reference model corresponding to the code and progenitor (e.g. s20-pp-str-mb-norel is mapped from s20-ref-1D; cf. Table 1 and Sect. 2.3) – hence the collapse and bounce are always simulated without any transport simplifications. We note, however, that focussing here only on the post-bounce evolution does not mean that any of the considered modeling variations is insignificant when applied during the collapse phase. In fact, some variations may even lead to more dramatic consequences than observed here when used during the collapse (e.g. Arnett, 1977; Bruenn, 1985; Lentz et al., 2012b).
We employ the equation of state (EOS) “SFHo” of Steiner et al. (2013). Since this EOS is constructed assuming nuclear statistical equilibrium (NSE), ideally we should replace it by a composition-dependent EOS linked to a nuclear reaction network once the temperature drops below GK. However, for the present study we avoid this additional level of complexity and employ the SFHo-table everywhere. Any neutrino interactions with light nuclear clusters (i.e. elements with mass numbers except particles) appearing in the SFHo EOS are ignored. The original SFHo table only allows electron fractions up to . This turned out to be insufficient for a subset of models in which transient, proton-rich bubbles arise in the gain region. To this end, we extended the original SFHo table in the sub-nuclear domain towards more proton-rich conditions using a 23-species NSE solver.
For the Alcar simulations, the radial grid remains fixed during the simulations. The standard resolution has radial zones with a width of m below km and growing by a constant factor per cell up to km at km and finally increasing by per cell up to the outer grid boundary at cm. In the Vertex simulations the radial grid at the time of bounce consists of zones, the distribution of which results from a quasi-Lagrangian treatment of the collapse, while during the post-bounce evolution the radial grid remains fixed (i.e. Eulerian) but is manually refined multiple times around the neutrinosphere by using a density-based criterion. Typically, a final number of radial zones of is reached at the end of the simulations. The angular grids are always uniform in , with the standard number of zones being and for Alcar and Vertex, respectively. The neutrino energy space is discretized using 15 energy groups distributed nearly logarithmically between 0 and MeV (Alcar) or MeV (Vertex).
With a spherical polar coordinate mesh the Courant time step would become prohibitively small near the coordinate center. We mitigate this problem by using a spherically symmetric core up to a certain radius. In Alcar, this radius decreases linearly in time from 10 km at 20 ms post bounce down to 7 km at about 0.5 s post bounce, whereafter it remains constant. In Vertex, the radius of the 1D core is 1.6 km at all times.
Our reference models include the following set of neutrino reactions (and corresponding reverse reactions):
Electron captures by heavy nuclei as in Bruenn (1985).
Axial response corrections to the neutrino-nucleon scattering cross section due to many-body effects as in (Horowitz et al., 2017, using the fit formula for ).
However, the impact of these opacity upgrades on the post-bounce evolution has been discussed before, so the main reasons for including them here are, first, to provide a reference for comparison with the considered modeling simplifications, and second, to enhance for some models the explodability in order to test modeling simplifications in an extended range of conditions (see Sect. 2.3 for the investigated list of models).
For the implementation of the rates in Alcar we followed as closely as possible that of Vertex, which is described in the Appendix of Rampp & Janka (2002). However, all neutrino rates (as well as all other physics modules) in Alcar have been coded from scratch, i.e. no routines have been copied from Vertex. The numerical treatment of the source terms in Alcar differs from that of Vertex in that it, first, allows for three instead of only one flux-vector component, and second, is both explicit and implicit in time depending on the type of interaction and dynamic conditions, whereas it is fully implicit in Vertex. We provide more details on the implementation of the source terms in Alcar in Appendix A.
We note that the set of neutrino interactions chosen here is not as advanced as the one usually employed in Vertex (e.g. Buras et al., 2006). However, we chose this rather basic set to facilitate the comparison between the two codes and to allow other groups to compare with our results on the basis of a widely available repository of input physics.
2.3 Investigated models
The axisymmetric reference models run with Alcar are s20-ref and s9-ref for the s20 and s9 progenitors, respectively, and evolve the fully multidimensional (i.e. not RbR-constrained) transport equations, Eqs. (3), using the aforementioned neutrino interactions (except strangeness and many-body corrections). Models with “str” and “mb” in their names additionally include the strangeness and many-body corrections, respectively. Models with “rbr” in their names adopt the RbR+ approximation in Alcar exactly as described in Buras et al. (2006), namely by setting the lateral flux densities444We note that this also includes the diffusive component of the HLL fluxes through the cell interfaces, which can be non-zero even if the cell-centered fluxes vanish., (though is fulfilled anyway in our non-rotating models) and setting the lateral source term for the evolution of , , to
for trapping densities, g cm, and to zero for lower densities. In this way, lateral advection and compression of energy and radial flux remain included in the transport equations as well as lateral neutrino-pressure forces in the gas-momentum equations.
The remaining modeling variations are: Turning off inelastic neutrino-electron scattering (“nones”), ignoring velocity-dependent and gravitational redshift terms (“norel”) by setting and in Eqs. (3), and using the simplified description of (-annihilation and bremsstrahlung) pair processes (“pp”) as suggested in O’Connor (2015). The pair-process simplification consists of, first, ignoring all pair processes for electron-type neutrinos entirely, and second, assuming that respective annihilation partners for pair-annihilation of are in isotropic, local thermodynamic equilibrium (LTE) when computing the pair-process interaction kernel. The second assumption reduces the source terms for pair processes to be formally equivalent to source terms for emission/absorption processes (see O’Connor, 2015, for more details). All three aforementioned simplifications are particularly appealing for energy-dependent transport schemes, because each of them entails dropping numerically complicated energy-bin coupling terms. We include these simplifications here, because they have rarely been tested so far in multidimensional simulations555One exception is inelastic neutrino-electron scattering, which was found in Burrows et al. (2018) to have a similar impact as seen in this paper. For 1D models the impact of neglecting neutrino-electron scattering and velocity terms has been tested in Lentz et al. (2012b, a). However, those studies did not test the impact just on the post-bounce evolution (i.e. starting from the same post-bounce models), which is what we do here..
The list of all simulations is provided in Table 1. The two Alcar models using RbR+ with and without neutrino-electron scattering (s20-rbr and s20-rbr-nones, respectively) are compared with the two Vertex models s20VX and s20VX-nones that contain exactly the same corresponding input physics. The setup of the remaining models is motivated mainly to test the impact of each modeling simplification (represented by “rbr”, “nones”, “norel”, “pp”) on the eventual onset time of explosion, and to compare this impact with that of using the opacity improvements labeled by “str” and “mb”. For the s20 model666The reader might wonder why for some of the one-dimensional s20 models the corresponding 2D counterparts are missing. The reason simply is that these counterparts are too unlikely to explode given the existing results, and therefore were deemed less informative for this study., each modeling simplification is tested for at least two models with different conditions regarding the proximity to explosion; e.g., we turn off neutrino-electron scattering both for rather late (s20-rbr/-nones) and early (s20-rbr-pp/-nones) exploding models. Since the reference s20 model, s20-ref, does not lead to an explosion, we use in some cases its exploding RbR+ counterpart, s20-rbr, as reference for the comparison. Additional models (s20-pp-str, s20-pp-mb, and s20-pp-str-mb) are set up to test which and how many opacity variations are needed to push the reference s20 model to an early explosion.
In order to obtain a rough idea about the influence of stochasticity, we perform for some models additional simulations that differ only in the pattern (but not amplitude) of initial random density perturbations. These simulations are labeled by numbers at the end of the model names. When discussing these models below we will suppress these numbers if the point of concern holds for all simulations of the given model independently of the initial perturbation pattern. Finally, several simulations with different numbers of grid cells are conducted in order to test the resolution dependence.
All simulations are stopped either once shock expansion sets in or after 1 s (0.9 s) of post-bounce evolution in Alcar (Vertex), except for individual models that are followed slightly longer because of optimistic runaway conditions.
3 Results: 1D models
We begin by comparing spherically symmetric models. Although rather unspectacular, 1D models offer the most straightforward and computationally economic way to identify the basic impact of modeling variations, and they may reveal features that could remain obscured by stochasticity in multidimensional models. Since for non-exploding 1D models the progenitor dependence is less relevant, we only discuss the s20 models in this section but not the s9 models.
3.1 Comparison with Vertex
The phases of collapse and early post-bounce evolution are simulated for all 1D and 2D models with the same input physics, thus we only need to compare these phases for models s20-ref-1D and s20VX-1D. The main features of the collapse are displayed in Fig. 1, which shows the electron-fraction, , lepton fraction, , and entropy per baryon, , in the stellar center as functions of the time-evolving density, , at the same location, as well as profiles of and along the enclosed mass coordinate . The core deleptonization and heating, aided by neutrino-electron scattering (e.g. Bruenn, 1985; Thompson et al., 2003), proceeds up to the trapping density of g cm, above which and remain essentially constant. The central values as well as the profiles agree well for both codes. The bounce777By definition, the bounce happens at the time when the shock entropy reaches 3 baryon. commences ms (298 ms) after initialization for the Alcar (Vertex) model at a mass coordinate of (). In Fig. 2 we show profiles of various quantities as functions of the enclosed mass (for times at and short after bounce) and radius (for later times). As can be seen in Fig. 2, the profiles of hydrodynamic quantities agree almost perfectly between Alcar and Vertex within the first tens of milliseconds after bounce. At this early stage the only noteworthy differences between both codes are observed in the neutrino-related quantities, e.g. in the luminosity,
and energy-averaged Eddington factor,
a few milliseconds after bounce right when the burst of electron neutrinos is released by the outward traveling shock wave. The aforementioned differences are related to what is seen in the small inset of the top left panel in Fig. 4, namely that the burst in Alcar has % higher peak luminosity (erg s vs. erg s) and contains % more energy (erg vs. erg as integrated between ms and ms) than the burst in Vertex. This difference goes hand in hand with a transiently lower deleptonization in the Vertex run (see Fig. 2, left column, third panel at ms). The reason remains unknown; it might be related to the approximate two-moment closure for the transport employed by Alcar. The impact of the more energetic burst on the subsequent evolution is difficult to quantify, but it must be relatively small given that the agreement between both codes remains very good at later times. We also point out that Vertex exhibits a somewhat stronger diffusive broadening of the outward propagating burst compared to Alcar. This may be connected to a finer radial grid, smaller time step, and higher-order spatial reconstruction scheme used for integrating the transport equations in Alcar. This difference in the translatory behavior of the burst at large distances has no relevance for the model evolution.
Before comparing the post-bounce dynamics below the shock we take a look at the mass accretion rates measured at km in Fig. 3. The basically perfect match between Alcar and Vertex provides confidence that the numerical treatment of gravity and the tabulated EOS is handled consistently in both codes and, equally important, ensures that any differences of dynamics seen below the shock are not the result of different conditions above the shock.
In Fig. 4 we show important global quantities as functions of time for the four Alcar and Vertex models with and without scattering. The trajectory of the shock888Since the shock is not a perfect discontinuity but is distributed over several radial zones, we compute the shock radius as the arithmetic average of the radii bracketing the extended, numerical shock., , given in the top right panel of Fig. 4, reveals that prompt shock expansion occurs for about 70-80 ms, whereafter the shock stalls and recedes. Not surprisingly, in the present 1D models the shock retraction takes place smoothly and monotonically, except for a bump at s that is related to the drop in the mass accretion rate (cf. Fig. 3) due to the infalling Si/Si-O interface.
The neutrino luminosities, , mean energies999We note that the mean energies computed as in Eq. (8) are at km essentially identical to the mean energies of the neutrino flux (obtained through replacing in Eq. (8) by and by ) because is very close to unity at this large radius.
(both plotted in Fig. 4 at km in the lab-frame as seen by an observer at infinity), the root-mean-squared (rms) energies,
(plotted in Fig. 4 at the gain radius as seen by an observer locally comoving with the fluid), neutrino heating rates in the gain layer,
(with luminosities defined as in Fig. 4), as well as the characteristic timescales of neutrino heating and fluid advection through the gain layer (with mass and sum of kinetic, internal, and gravitational energy of ),
respectively, agree for most of the time to within a few per cent between both codes101010We adopt the same definition of the aforementioned quantities as Summa et al. (2016), except for possibly different locations and reference frames of measurement and the fact that the rms-energies here are based on the energy distribution, , instead of the number distribution, , such that the global heating rate can be approximately written as , where is the gain radius. Since only electron-type neutrinos contribute significantly to gain-layer heating, we omit plots and discussions of for neutrinos..
Within the small level of disagreement we notice the tendency of both Alcar models with and without neutrino-electron scattering to produce slightly higher neutrino luminosities, , in both electron-neutrino species. This enhancement of energy-loss rates in Alcar may be the reason why we observe a slightly smaller radius, , and higher temperature, , of the PNS surface (defined by the location where g cm; cf. Fig. 4), higher mean energies, , for all neutrino species, a less extended shock radius, , and gain radius, , as well as a higher mass-specific neutrino heating rate, , for the Alcar models compared to the Vertex models. Remarkably, however, the ratio of the characteristic timescales, (shown in the bottom left panel in Fig. 4), which is a more meaningful measure of the proximity to explosion than the aforementioned quantities (e.g. Janka, 2001), does not exhibit a clear trend in either direction. This means that the impact of a more compact configuration in Alcar, which in the first place is detrimental to an explosive runaway because the gain layer and shock sit deeper in the gravitational well, is approximately compensated by higher luminosities and therefore more powerful neutrino heating.
3.2 Electron scattering
Figure 4 also shows the corresponding 1D results obtained when switching off neutrino-electron scattering at about ms after bounce for both Alcar and Vertex. In the present 1D models the impact on many quantities is only on the percent level. Nevertheless, the less efficient energy deposition without neutrino-electron scattering (cf. in Fig. 4) is sufficient to reduce the shock radius, , by a few km and to lead to a more sizable reduction of the advection timescale, , and of the timescale ratio, , by . The most notable (but for the dynamics rather irrelevant) difference concerning the neutrino emission appears for the mean energies, , of the heavy-lepton neutrinos, which are several MeV higher without the efficient down-scattering process of neutrino-electron scattering (Fig. 4, left column, second panel from top).
3.3 Pair processes
We start by considering model s20-pp-1D that incorporates the simplified pair processes treatment, i.e. for electron-type neutrinos all pair processes are neglected and for pair-annihilation of neutrinos the corresponding annihilation partners are assumed to be in isotropic LTE. The top left plot in Fig. 5 shows that the luminosities of heavy-lepton neutrinos are reduced compared to the reference case by during the first ms of post-bounce evolution. The main reason for this reduction is most likely the approximate assumption that pair-annihilation targets are isotropically distributed in momentum space, whereas they actually become more and more forward peaked with increasing radius. This boosts the pair-annihilation rates while leaving the rates of the inverse (i.e. pair-production) reactions unchanged. We note that an impact of similar size has been found also by O’Connor (2015) for 1D models with electron-positron annihilation but without bremsstrahlung. Keeping in mind that four times the individual luminosity enters the total energy loss rate of the PNS, this reduction of neutron-star cooling during ms probably explains the observed increase of the neutron-star radius by km and of the shock radius by km compared to the reference model.
The luminosities and spectral properties of emitted electron-type neutrinos remain fairly unaffected by the pair-process simplification. This helps understanding the similarly weak sensitivity of the specific heating rate, which can approximately be written as
(where ), and the heating timescale, which roughly scales like
(where we assumed the specific total energy of the gain region to be approximately proportional to the gravitational energy at the gain radius). The heating rate, heating efficiency and the advection timescale, on the other hand, show stronger deviations from the reference model, mainly because they are directly proportional to the mass in the gain layer, which itself is enhanced by compared to the reference model. As a result, , which is most relevant for the explosion behavior, is enhanced by a comparable amount during almost the entire evolution.
We note that this result is not in tension with previous studies (and with the results for models s20-str-1D and s20-mb-1D discussed below) that report more favorable runaway conditions for cases in which the PNS contraction proceeds faster (e.g. Marek et al., 2009; O’Connor & Couch, 2018; Melson et al., 2015b; Bollig et al., 2017). This is because in those studies the faster PNS contraction is the result of changing physics ingredients different from the ones varied here, e.g, using a softer nuclear equation of state, general relativity instead of Newtonian gravity, or reduced nucleon-scattering opacities for all neutrinos. In those cases the accelerated PNS contraction came along with more favorable neutrino emission properties (i.e. higher values of ) that overcompensated for the stronger gravitational binding. In the present case, on the other hand, the deceleration of PNS contraction due to artificially amplified annihilation of pairs is not accompanied by a large enough reduction of to cause more pessimistic runaway conditions.
3.4 Strangeness and many-body corrections
The two 1D models, s20-str-1D and s20-mb-1D, both incorporate corrections that effectively reduce the opacities for neutrino-nucleon scattering, in particular in the crucial neutrino-decoupling region, by a few percent. As a result, we observe an enhancement of the neutrino luminosities for all three species by a similar amount, which is well in agreement with previous studies that incorporated these corrections (e.g. Melson et al., 2015b; O’Connor et al., 2017). In the present one-dimensional case, this leads to smaller shock- and PNS-radii throughout, but slightly higher mean- and rms-energies of neutrinos and therefore, in combination with the increased luminosities, to slightly stronger neutrino heating, shorter and higher , thus slightly more favorable conditions for shock revival.
3.5 Velocity-dependent and gravitational redshift terms
For model s20-norel-1D, in which frame-dependent effects such as advection and Doppler-shift as well as gravitational redshift are ignored after ms, one immediately recognizes the significantly smaller heating rate as well as specific heating rate, reduced heating efficiency, longer heating timescale, and reduced timescale ratio compared to the reference case. At first glance, this seems counterintuitive in view of the fact that according to Fig. 5 the (lab-frame) luminosities agree well and the (lab-frame) mean energies are even higher compared to the reference model. However, both quantities are not sufficiently representative of the heating conditions in the gain layer, as is revealed by Fig. 6 that depicts radial profiles of the luminosities and rms-energies measured in the comoving frame111111For model s20-norel-1D the distinction between the comoving frame and lab frame is meaningless and the evolved neutrino moments are used for plots of both comoving- and lab-frame quantities. at a representative post-bounce time for both 1D and 2D models: For model s20-norel-1D the neutrino fluxes as seen in the fluid frame are almost smaller in the gain region than for the reference model. The good agreement between the luminosities of s20-norel-1D with the local lab-frame luminosities of the reference model (dashed lines in Fig. 6) suggests that this difference stems from Doppler boosting in the infalling material, which is ignored in model s20-norel-1D. One might still wonder why the rms-energies, (cp. Figs. 5 and 6) exhibit hardly any reduction in model s20-norel-1D compared to the reference model. The most likely explanation is that the frame correction is in relative terms smaller for the rms-energies (transforming as for forward peaked radiation) than for the luminosities (transforming as for forward peaked radiation), such that gravitational redshift may approximately compensate for the effect of Doppler blueshift for the rms-energies.
4 Results: 2D models
After having obtained an idea of the level of agreement between the two codes Alcar and Vertex and the impact of our modeling variations for the spherically symmetric case, we now investigate how these dependencies translate to the 2D axisymmetric case. Additionally, we will examine the impact of using the RbR+ approximation for the Alcar code, comment on the level of stochasticity and numerical convergence, and contrast our results for the s20 models with others found in the literature.
4.1 Basic features
Before comparing the models in detail, we first summarize some features common to most models. Figure 7 depicts snapshots of the isotropic-equivalent lab-frame luminosity, , and the specific entropy taken at three post-bounce times for two Alcar models and one Vertex model, and Figs. 8, 10, and 11 provide a summary of global properties for most investigated models. By switching to 2D, we allow the fluid to develop fluid instabilities such as PNS convection, SASI, and post-shock convection driven by neutrino heating.
Proto-neutron star convection is the consequence of unstable gradients of the lepton number and specific entropy developing below the PNS surface as a result of neutrino emission. Leptons trapped in the dense core of the PNS are shuffled into the overlying layers, where they provide additional pressure support and thereby slow down the PNS contraction. Correspondingly, we observe larger PNS radii, , and lower temperatures, , as well as reduced energies of emitted neutrinos, and , in 2D compared to 1D at equal times (Buras et al., 2006; Dessart et al., 2006). The luminosities, , of electron-type neutrinos remain almost unchanged, while those of the heavy-lepton neutrinos are enhanced by several tens of percent. All our models (although with some differences for those models neglecting velocity terms in the transport, see Sec. 4.4), reproduce these effects connected to PNS convection in good agreement with previous studies (e.g. Buras et al., 2006; Dessart et al., 2006; Müller et al., 2012a; Bruenn et al., 2016; Radice et al., 2017).
While PNS convection already starts at about ms after bounce, the gain layer remains fairly spherically symmetric for a much longer time, namely until about s (s) for the s20 (s9) progenitor models, as can be recognized in Figs. 8 and 10 (11) for the s20 (s9) models by the correspondingly late rise of the lateral (i.e. carried by non-radial motions) kinetic energy integrated over the gain layer, . This means that until these times the conditions for efficient growth are met neither for SASI nor for post-shock convection. Since the hydrodynamic instabilities are triggered by random perturbations, the exact time of significant departure from spherical symmetry slightly varies between the models. It may be worth pointing out that in tests with Alcar using the s20 model we experienced quite some sensitivity of this transition time to the numerical treatment: Using a low (i.e. linear) order for spatial reconstruction of the hydro variables or employing the more diffusive HLL solver led to a significantly earlier onset of non-radial flow activity.
4.1.1 s20 progenitor models
Since SASI growth is favored by short advection timescales (Blondin & Mezzacappa, 2006; Foglizzo et al., 2007; Scheck et al., 2008), , efficient growth sets in once the post-shock configuration has become sufficiently compact. Around ms after bounce the exponential growth of SASI modes with periods comparable to can be seen in the dipole component of the shock surface (cf. bottom left panel of Figs. 12), which is obtained from a decomposition of the latter into Legendre polynomials (see, e.g., Summa et al., 2016, for the exact definition of the multipole coefficients). The snapshots in the top row of Fig. 7 are taken at ms right around the time when the transition from the well-ordered, linear phase to the turbulent, non-linear phase of the SASI takes place. After this transition the evolution of the post-shock layer is characterized by parasitic instabilities (Rayleigh-Taylor and Kelvin-Helmholtz) whose turbulent mass motions tap energy from the SASI. Post-shock convection tends to be suppressed for non-exploding s20 runs because of the retreating shock radius and correspondingly short advection timescales, which act against the development of a sufficiently large negative entropy gradient and correspondingly short growth timescales for convection. This statement is backed by the circumstance that the -parameter characterizing the growth conditions of post-shock convection (see, e.g., Foglizzo et al., 2006; Summa et al., 2016 for the definition of this quantity and details regarding its interpretation and computation) remains below the critical value of (bottom right panel of Figs. 8 and 10). Nevertheless, since the critical condition for convection, , holds strictly only for the linear regime, the existence of secondary convective instability associated with highly non-linear SASI activity is not excluded.
4.1.2 s9 progenitor models
For the models using the s9 progenitor, which are characterized by considerably smaller mass accretion rates with respect to the s20 models (see Fig. 3) and therefore a less compact shock surface, the situation is reversed in that the dominant fluid instability is not SASI but post-shock convection. Correspondingly, we observe a transition of the -parameter above the critical value of at about s (cf. Fig. 11), whereafter remains during the entire simulation. This strong progenitor dependence of instability regimes, with massive (low-mass) progenitors favoring SASI (convection) has been recognized before, e.g. in Müller et al. (2012b); Fernández et al. (2014).
4.2 Comparison with Vertex
We start by comparing the two Alcar models that incorporate the RbR+ approximation with (s20-rbr) and without (s20-rbr-nones) neutrino-electron scattering to the corresponding Vertex runs (s20VX and s20VX-nones, respectively).
First of all, as visible in Fig. 8 and Table 1, the models without neutrino-electron scattering do not explode (at least not until the end of each simulation), while the ones including neutrino-electron scattering do explode, but rather late (i.e. several hundred milliseconds after the infall of the Si/Si-O interface at ms). The agreement of both codes in clearly showing the impact of a relatively small variation of neutrino-interaction rates (in the present case due to neutrino-electron scattering) is thus already encouraging. Moreover, for both codes the exploding models are characterized by a significant scatter in explosion times (cf. Table 1): The three Alcar models with the same input physics but different random initial perturbation patterns cover a large range of explosion times of s, while the two corresponding Vertex runs explode also within this time interval, namely at s and s.
The overview of important global properties as functions of time in Fig 8 reveals that the differences between both codes concerning the neutrino luminosities and energies, as well as the PNS surface temperature and radius remain on the few-percent level, as found in 1D. This suggests that PNS convection is operating in both codes consistently regarding its impact on PNS contraction and neutrino emission. Nevertheless, a closer look unfolds that the and luminosities now show a small enhancement in Vertex relative to Alcar for s that was not seen in 1D (see Fig. 4) and may therefore be connected to some discrepancy in the PNS convection. Indeed, an inspection of the radial profiles of the lateral kinetic energy densities and the neutrino luminosities in Fig. 9 reveals that the PNS convection in Alcar proceeds somewhat differently with higher kinetic energies. Unfortunately, despite a dedicated analysis we were unable to track down the exact reason for this discrepancy and its detailed consequences. We can, however, already clearly assess that this enhancement in convective energy is not an artifact of the RbR+ approximation, because model s20-ref1, which does not employ RbR+, shows a similar behavior (cf. Fig. 9 and a more detailed discussion of the RbR+ approximation in Sec. 4.5). It may be that the stronger PNS convection in Alcar is related to a subtle difference in the PNS structure, which is systematically more compact (km difference in ) than in Vertex models in 1D (Fig. 8) as well as 2D. This might allow PNS convection to be enhanced in Alcar by tapping the higher gravitational binding energy of the more compact PNS.
We do not deem the disagreement in the PNS convection zone to be overly significant concerning the explosion dynamics, because the shock radius, the neutrino heating rates, characteristic timescales, as well as most remaining properties are in good agreement apart from stochastic fluctuations.
Assessing in detail the heating conditions in the gain layer and the shock evolution, and ultimately finding the reason why a simulation at a certain time exhibits shock expansion or not, is difficult given the complicated temporal behavior of heating-related diagnostic quantities such as , , , and . At later times, s, we observe, consistently for both codes, that the exploding models in comparison to the non-exploding models exhibit temporal variations on longer timescales and with larger amplitudes. For instance, , occasionally travels up and down between on timescales of hundreds of milliseconds for the exploding models, while for the non-exploding models it remains rather close to with relatively short and low-amplitude oscillations. The larger amplitudes and longer timescales of temporal variations suggest that the exploding models hover in a state very close to criticality before they ultimately explode. Retaining in such a state means that small perturbations may have a large dynamic impact. A natural suggestion from this is that the temporal pattern observed for the heating-related properties is shaped more by stochasticity than systematics. This conjecture is supported when comparing models that were initialized with a different pattern of density perturbations (cp. thick and thin lines of same color in Fig. 8): , for instance, varies over a comparable range of values but with a considerably different temporal behavior for different initial perturbations. We infer that the time of shock runaway must be similarly affected by stochasticity, which is confirmed by the large scatter of runaway times seen for both codes. Therefore, even though the explosion times of individual simulations may differ by several hundred milliseconds, we conclude that the overall agreement between both codes is good and no worrisome differences concerning the explosion behavior are found.
Before moving on to the next topic, we comment on a rather peculiar difference between Alcar and Vertex, namely that in Alcar the evolved neutrino energies and fluxes appear to have rather fine spatio-temporal structures while the neutrino quantities are smoother in Vertex (see, e.g., the color maps of the isotropic-equivalent luminosities in Fig. 7). We suspect two possible reasons: First, the approximate, non-linear closure in Alcar, which allows self-interaction and non-linear effects (such as shocks) of the radiation moments even in optically thin regions (e.g. Pons et al., 2000), and second, the use of an explicit time integrator in Alcar as opposed to an implicit one in Vertex that tends to broaden sharp features. At this point it is thus not clear up to which level the observed small-scale features in Alcar are physical or not; however, the dynamical consequences seem to be marginal given the otherwise good agreement.
4.3 Electron scattering, pair processes, strangeness and many-body corrections
As already pointed out in the previous section the relatively small effects associated with neutrino-electron scattering, which has often been neglected in a number of previous studies assuming that it would only be relevant during collapse, turn out to be a crucial ingredient for initiating shock runaway for the s20 model (in qualitative agreement with Burrows et al., 2018). For the s9 progenitor models (see Table 1 and Fig. 11 for an overview of global properties as functions of time) we see the same strong sensitivity, i.e. models exploding around s (s9-ref, s9-rbr) are turned into non-exploding models (s9-nones, s9-rbr-nones) when ignoring neutrino-electron scattering. These results demonstrate once more how sensitive the onset of shock runaway can be to variations in the neutrino interaction rates.
In view of the substantial impact of neutrino-electron scattering it seems not too astonishing that we see a similarly strong impact of other variations of the neutrino rates considered in this study, both for the s20 (see Fig. 10) and s9 (see Fig. 11) models: The simplified treatment of pair processes (used in all models labeled by “pp” in their names), as well as the scattering-opacity reducing strangeness and many-body corrections (including “str” and “mb” in their names, respectively) are all conducive to explosions, i.e. each modification comes along with a visible boost of the timescale ratio, , as well as (in the case of a successful explosion) a shift towards earlier times of shock expansion. This shift is naturally larger in the s20 models, because the reference s20 model (s20-ref) lacks an explosion while the reference s9 model (s9-ref) already explodes around s. Combining all three of these variations turns a robustly (in the sense of being independent of stochasticity) non-exploding model (s20-ref) into a robustly exploding model (s20-pp-str-mb), with explosion times between 0.38-0.50 s, cf. Table 1.
Crudely judging from the explosion times (cf. Table 1), the impact of the pair-process simplification is roughly comparable to that of the combination of the two physically motivated opacity corrections (“str” and “mb”; cp. models s9-ref, s9-str-mb, and s9-pp) and slightly stronger than that of neutrino-electron scattering (cp. models s20-rbr-pp, s20-rbr, and s20-rbr-pp-nones).
Finally, we point out an interesting feature: The impact of each microphysics variation on the explosion behavior is qualitatively correctly predicted in 1D by the timescale ratio, , in the sense that all modifications, which lead to larger values of in 1D, trigger an earlier onset of explosion, or an explosion at all, in 2D (see for the s20 models Figs. 5 and Sec. 3, as well as for the s9 models the dotted lines in the panel for in Fig. 11 showing the behavior of the corresponding 1D models). In contrast, the shock trajectory of 1D models is not a good indicator of more/less favorable runaway conditions: Both strangeness and many-body corrections reduce the shock radius in 1D, while the pair-process simplification increases the latter (cf. Fig. 5), but all three of these physics variations help initiating a shock runaway in 2D.
4.4 Velocity-dependent and gravitational redshift terms
Concerning the prospects of shock runaway, the purely negative impact of neglecting velocity-dependent as well as gravitational redshift terms in the transport seen in 1D seems to be partially alleviated in 2D, albeit in a model-dependent fashion: An early exploding s20 model (s20-pp-str-mb) explodes even earlier with these modifications (s20-pp-str-mb-norel), a late exploding s20 model (s20-rbr) then lacks an explosion (s20-rbr-norel), and the reference s9 model (s9-ref) seems almost unaffected (s9-norel).
These results are most likely connected to the reduced convective activity in the PNS, which can be inferred from the smaller values of absolute velocities, , in the PNS convection zone displayed for a representative time in Fig. 6. The consequence is a PNS that contracts faster, is hotter, and therefore emits neutrinos with higher rms-energies (see , , and in Figs. 8, 10, and 11). For a less massive PNS these explosion facilitating consequences seem to be able to (over-)compensate for the missing Doppler blueshift in the infalling material, because the infall velocities are lower compared to a more massive PNS (see profiles of the radial velocity, , in Fig. 6). This dependence on the PNS mass at the time of close proximity to criticality might explain why the runaway conditions become more (less) optimistic for the originally early (late) exploding model when switching off the considered terms. In the models with the less massive s9 progenitor, all velocities, including those in the PNS convection zone, are lower to begin with, which might explain the almost vanishing net impact on the shock trajectory in this case.
The actual reason for the less efficient PNS convection when switching off the considered terms is not determined easily. We suspect the following: When neglecting velocity-dependent terms in the moment equations, neutrinos (i.e. their energy and lepton number) are not advected with the flow as they should, being deep inside the PNS. Instead, neutrinos carried in a bubble that is about to rise in response to an unstable stratification effectively leak out of the bubble and are left behind. Since near the bottom of the PNS convection zone the neutrino lepton number can be a sizable fraction of , these neutrino losses effectively reduce the lepton number of the bubble and therefore the tendency of the latter to rise further.
4.5 Ray-by-ray-plus approximation
In this section we compare the impact of using the RbR+ approximation on the PNS convection and neutrino emission as well as on the explosion dynamics.
4.5.1 PNS convection and neutrino emission
We start by addressing the question if using the RbR+ approximation might have a significant impact on the PNS convection and, in turn, on the neutrino emission from the deleptonizing PNS. Figure 9 compares the angle-averaged, radial profiles of (comoving-frame) luminosities, for all evolved neutrinos species, , as well as the radial profiles of the convective lateral kinetic energies, , at two representative times for models s20-ref and s20-rbr. The convectively unstable region is characterized by large values of . The neutrino luminosities – all the way from inside the PNS up to the saturation (i.e. gain) radius – as well as the convective energies are very similar for both models (apart from small temporal fluctuations) even at late post-bounce times. Moreover, as can be seen in Figs. 8 and 10 showing functions of time for the s20 models and Fig. 11 for the s9 models, also the mean and rms-energies of emitted neutrinos are barely affected. These results suggest that RbR+ has a negligible impact on PNS convection and on the angle-averaged neutrino emission, at least much less than, e.g., neglecting velocity-dependent terms in the transport (cf. Sec. 4.4). This result is not too astonishing considering that neutrinos deep inside the PNS are strongly coupled to the medium such that, first, the lateral advection fluxes included in the RbR+ approximation strongly dominate the lateral diffusion fluxes, and second, lateral neutrino-momentum transfer to the fluid is well described by using pressure derivatives associated with assumed isotropic neutrino distributions, as applied in the RbR+ treatment (Buras et al., 2006).
4.5.2 Impact on explosion behavior of s20 models
While the reference model, s20-ref, is not showing any indication of shock runaway during s of post-bounce evolution, the corresponding RbR+ model, s20-rbr (which was compared against the corresponding Vertex model in Sec. 4.2), does explode, though rather late and with a substantial scatter in the time of shock runaway. Likewise, switching to the simplified pair-process treatment, which itself promotes explodability, yields late explosions without (model s20-pp) and early explosions with (s20-rbr-pp) the RbR+ treatment. In our set of s20 models, the net effect of RbR+ regarding the explodability is comparable to using the pair-process simplification together with either the strangeness (s20-pp-str) or the many-body correction (s20-pp-mb) but smaller than using all three of those together (s20-pp-str-mb).
In order to understand what pushes the RbR+ models closer to criticality, we consider Fig. 12, which shows for several models the excess/shortage of neutrino fluxes, rms-energies, temperatures, and mass fluxes close to the poles (i.e. averaged over the surfaces of two opposite cones defined by and ) with respect to their averages over the full sphere, as well as the normalized Legendre coefficients for dipole and quadrupole deformation of the shock surface. All exploding and non-exploding s20 models with RbR+ (including the ones not plotted in Fig. 12) show a sustained polar enhancement of temperatures, mass-fluxes (albeit with large temporal fluctuations) and electron-type neutrino fluxes. At the same time, the RbR+ models come along with higher shock-oscillation amplitudes and therefore more optimistic runaway conditions.
Better insight can be gained from Fig. 13, where we compare for models s20-ref1 and s20-rbr1 color maps of several quantities averaged in time between s and s, namely the lateral velocity, , radial velocity, , and angular variations of the radial mass-flux density, , temperature, , and radial, energy-integrated flux density of electron-type neutrinos, . The angular variations are computed for each time average of quantity as
where denotes the angle average of at a given radius. For both models, the color maps of , and carry the imprint of a large-scale, low-order flow pattern, which in the present case is mainly the result of quasi-periodic SASI sloshing: At large radii, km, post-shock material preferentially expands and moves away from the symmetry axis while falling back down towards the PNS. This causes large amounts of gas to arrive at the cooling region in the vicinity of the equator (cf. near km in Fig. 13), while the sloshing motion is further carried on by material moving towards the correspondingly other hemisphere at radii km close to the PNS surface (e.g. for ). This portion of the flow converges near the poles and helps driving the next half-cycle of approximately cigar-shaped shock expansion. Importantly, the converging flow also results in a net accumulation of matter near the poles. As can be seen in the panels for by the red-colored polar regions around radii of km, this surplus of matter leads to (in absolute terms) increased mass-flux densities down onto poles of the PNS surface (the latter being roughly coincident with the largest purple circle denoting the g cm surface of the time-averaged configuration). These accretion hot spots produce enhanced temperatures at and high radial neutrino fluxes streaming away from the poles, consistent with the results of Fig. 12.
While the basic flow pattern and the polar temperature- and neutrino-flux enhancements are qualitatively similar in both models, the RbR+ model exhibits important differences compared to the unconstrained model: Non-radial motions are faster, polar expansion flows encountered during each SASI cycle are more powerful and even lead to positive time averages of radial velocities (red regions for km in the panel in Fig. 13), and the polar accretion hot spots are characterized by higher temperatures and stronger neutrino emission. Moreover, neutrino radiation clearly remains more collimated above the polar hot spots than in the model without RbR+, where the radial fluxes relative to their angle averages quickly decrease with increasing radius, probably as a result of high lateral neutrino fluxes pointing away from the hot spots.
Based on the results described above it seems reasonable to suspect that a positive-feedback mechanism may be at play between the originally purely fluid-dynamical advective-acoustic/advective cycle (e.g. Foglizzo, 2002; Blondin & Mezzacappa, 2006; Scheck et al., 2008; Guilet & Foglizzo, 2012) of the SASI and accretion-induced neutrino heating: Linear sloshing modes that lead to accretion hot spots at the PNS poles may get amplified due to absorption of neutrinos originating from these hot spots while, in turn, the neutrino emission rates at these hot spots may get boosted with growing amplitude of the sloshing modes. The efficiency of this feedback mechanism is higher with RbR+ than without, apparently because neutrino fluxes stemming from the hot spots remain more collimated when reaching the gain layer (see right panels in Fig. 13), causing more energy to be pumped into the fastest, near-axis material that carries most of the kinetic energy.
If our picture is correct, we may additionally speculate that the rather significant impact of RbR+ in the present SASI-dominated models is fostered by the following circumstances: First, sloshing modes generate strong downflows always at the same location, namely near the poles, which facilitates the development of hot spots and creates locally enhanced neutrino fluxes that are always pointed towards the optimal direction for dynamical feedback. Second, sloshing modes cause downflows to be rather well synchronized with expansion flows, i.e. both take place on the same characteristic timescales (namely the advection timescale ) and downflows are always quasi-periodically succeeded by expansion flows.
4.5.3 Impact on explosion behavior of s9 models and comparison with s20 models
We now consider the convection-dominated s9 models; see Figs. 11 and 14 for the corresponding time-dependent properties of selected models, as well as Fig. 15 for color maps of time-averaged quantities for two s9 models without and with RbR+. Here, using RbR+ advances the onset time of shock runaway only by s in model s9-rbr compared to s9-ref (cp. Table 1), which is a rather small effect considering that neglecting neutrino-electron scattering prohibits the runaway entirely, both with (s9-rbr-nones) and without (s9-nones) using RbR+. The impact of RbR+ is thus much less extreme than in the s20 models, where switching off RbR+ delays the explosion by a much longer time than switching off neutrino-electron scattering (cp. models s20-rbr-pp, s20-pp, and s20-rbr-pp-nones, respectively).
The reduced impact of RbR+ for the s9 models becomes particularly obvious when comparing the two longest evolving, non-exploding models s9-nones and s9-rbr-nones (cf. Figs. 11, 14): The polar neutrino emission is only marginally enhanced by RbR+, and the dipole and quadrupole shock-deformation modes, the polar temperatures and mass accretion rates at the PNS surface, as well as essentially all global properties in Fig. 11 are nearly identical for both models.
What could be the reason for the much smaller impact of RbR+ in the s9 models compared to the s20 models? At least to some degree the answer must be connected to the fundamental difference between the two fluid instabilities respectively at play: The SASI is a global instability that triggers quasi-periodic, linear sloshing motions between the two hemispheres. Neutrino-driven convection, on the other hand, instigates buoyancy modes (i.e. high-entropy bubbles) on, at least initially, smaller spatial scales than SASI, corresponding to higher multipole orders (Foglizzo et al., 2006), and in a spatio-temporally much more stochastic fashion. It therefore seems plausible to suspect that large-scale asymmetries (associated with fluid modes of low multipole order) are typically smaller for a convection-dominated flow and, recalling the results for the s20 models, that a feedback mechanism between polar hot spots and neutrino heating is less likely to be realized.
However, this rather qualitative argument alone may be unsatisfactory without having considered in more detail the actual flow pattern resulting in the saturated state (cf. Figs. 14 and 15). For this to be properly understood we recall, first, that axisymmetric models exhibit a designated axis parallel to which bubbles tend to rise more readily than in any other direction, and second, that they are subject to an inverted turbulent energy cascade, i.e. energy is transported from small to large scales and not the other way around as in 3D (Kraichnan, 1967; Hanke et al., 2012; Couch, 2013). Consistent with these properties we observe that also the s9 models develop powerful low-order modes: The coefficients of dipolar () and quadrupolar () shock deformation modes in Fig. 14 reach values as high as or even higher than those of the s20 models (cp. Fig. 12). Moreover, the panels for , , and in Fig. 15 indicate that the time-averaged flow is dominated by a single large-scale convective eddy in each hemisphere: Hot bubbles rise preferentially near the poles, expand sideways while cooling down, re-enter the cooling layer near the equator, replace previously ascending bubbles at the poles, and finally start to rise again.
We caution the reader, however, that the circumstance that some qualitative features of the time-averaged flow appear to be similar in Figs. 15 and 13 does not imply similar dynamical behavior of both types of models, s9 and s20, respectively. On the contrary: In the convection-dominated s9 models the flow pattern is rather quasi-stationary and in each hemisphere largely decoupled from that in the correspondingly other hemisphere. In contrast, in the SASI-dominated s20 models the flow pattern is much more variable with time, because post-shock material quasi-periodically sloshes from one hemisphere to another.
Coming back to the question why the impact of RbR+ is reduced in the s9 models: The panels for and in Fig. 15 show that even though the s9 models do exhibit vigorous low-order (i.e. dipole and quadrupole) convective modes, the mass-accretion rate in this progenitor is so low that the convective flow does not accumulate enough matter near the poles of the PNS surface to lead to dynamically relevant emission hot spots.
Considering only two progenitor models here we can hardly make general statements concerning the impact of RbR+ in axisymmetry that hold for all convection-dominated and SASI-dominated models, led alone the cases where such a distinction is not possible. However, our results may be indicative of the tendency that a convection-induced flow pattern composed of two quasi-stationary eddies is generically not as efficient as SASI sloshing in producing – around the poles or at any other location on the PNS surface – a significant surplus of matter that could lead to a dynamically relevant emission asymmetry. Perhaps an additional systematic reason speaking against the possibility of a positive-feedback cycle in convection-dominated models may be the lack of synchronicity between down- and expansion flows, i.e. the circumstance that the rise time of convective bubbles (roughly given by ) is too long compared to the advection timescale, , for bubbles to experience efficient feedback from neutrino emission that is released by downflows on an advection timescale. In other words, the enhancement of polar neutrino emission triggered by some polar downflow declines too quickly for the next convective bubble to reach appreciable positive radial velocity. This argument could at least explain why for the two exploding models, s9-rbr and s9-ref, we actually do see an, although small, explosion-promoting impact of RbR+, because here is obviously closer to (cf. Fig. 11), at least shortly before shock runaway sets in.
4.6 Stochasticity and resolution dependence
We have already found (see, e.g., Sec. 4.2) that the impact of stochasticity and temporal fluctuations is low regarding the emission related properties but can be substantial for the gain-layer related properties including the onset time of explosion. In this section we want to briefly consider the results provided by our models that allow one to identify some systematic tendencies. Subsequently, we address the question how well the obtained global features, including the scatter of shock-runaway times, are converged with respect to the numerical grid.
Comparing the time evolution of (Figs. 8, 10 and 11) between various models reveals, first, that the amplitudes of temporal fluctuations are lower in the convection-dominated s9 models than in the SASI-dominated s20 models, and second, that at least for the s20 models these amplitudes grow with the proximity of each model to the runaway threshold. This tendency is consistent with the observed time interval, , within which the runaway times are dispersed for models differing only in the initial perturbation pattern (see Table 1): The scatter seems to be greater for the s20 models than for the s9 models, and for the s20 models the scatter is larger for more marginally exploding models (e.g. s20-rbr, s20-pp-str) than for more robustly exploding models (e.g. s20-rbr-pp). These tendencies are in qualitative agreement with Cardall & Budiardja (2015), who examined stochasticity using a large number of simplified models and found larger dispersion of explosion times for SASI-dominated models than for convection-dominated models (see also Kazeroni et al., 2017).
The strong dependence of on the progenitor and the input physics might be one reason why previous studies report quite diverse values of (e.g. O’Connor & Couch, 2018; Summa et al., 2016; Cardall & Budiardja, 2015; Takiwaki et al., 2014). An additional reason might simply be that many studies, particularly the ones including computationally expensive neutrino transport such as ours, are forced to rely on rather poor statistics because they can only afford a small number of simulations.
Coming now back to the second question concerning numerical convergence: For the reference s20 model, s20-ref, and its counterpart including the RbR+ approximation, s20-rbr, we repeated the simulations with both increased and decreased resolutions in both radial and angular directions (see models ending with “hires”, “lores”, “hi” and “lo” in Table 1), in some cases even multiple times with different initial perturbation patterns. We could not identify a systematic trend with varying the resolution, neither regarding the neutrino-emission properties nor the heating conditions in the gain layer nor the scatter in explosion times. In particular, for model s20-rbr the onset of explosion does not appear to be correlated with resolution in any direction: For all three angular resolutions with 80, 240, and 320 zones the scatter in the times of shock runaway remains comparably high, namely, s, 0.44 s, and 0.43 s, respectively. Since model s20-rbr is quite marginal concerning its tendency to explode and would therefore probably be quite sensitive to numerical resolution, this suggests (although does not prove) that, at least, features the runaway is sensitive to are numerically converged.
4.7 Previous axisymmetric simulations of the s20 progenitor
Here we want to collect some results concerning the explosion behavior of axisymmetric simulations performed by various other groups using the s20 progenitor. A detailed comparison is, however, out of the scope of this paper.
O’Connor & Couch (2018) used M1 transport in cylindrical coordinates, a similar neutrino setup as applied here for model s20-pp but additionally neglecting neutrino-electron scattering, and employing the LS220 EOS (Lattimer & Swesty, 1991). In their models shock expansion sets in around 700 ms post bounce, while our model s20-pp explodes later and would probably not explode at all without neutrino-electron scattering. However, the LS220 EOS used in O’Connor & Couch (2018) is slightly softer than the SFHo EOS used here and might lead to earlier explosion times.
Formally the same setup and a similar M1 scheme as in O’Connor & Couch (2018), although in spherical coordinates, was used by Skinner et al. (2016). We were unable to ascertain how Skinner et al. (2016) treated pair processes, e.g. if they ignored annihilation for electron-type neutrinos and used an LTE assumption for pair-annihilation partners as in O’Connor & Couch (2018) and as in our “pp”-models, or if they treated pair processes like we did in all other models. In the latter case, the lack of an explosion both with and without RbR+ would be consistent with our models, while in the former case it could mean that their models explode less readily than ours.
Kotake et al. (2018) employed an IDSA scheme in the RbR+ mode, the LS220 EOS, and a variety of different neutrino interactions, while energy-bin coupling reactions (neutrino-electron scattering as well as pair processes) are included up to 0th angular order of the interaction kernels. Their model G1 is most similar to our model s20-rbr-pp. While their model G1 does not explode until the simulation was stopped at 600 ms, our model explodes around 350 ms. However, the comparison is not conclusive because model G1 ignores weak magnetism and recoil corrections after Horowitz (2002).
Finally, Bruenn et al. (2016) and Summa et al. (2016), though using mutually different transport solvers, both employed RbR+ and an advanced (compared to the one used here) set of neutrino interactions and the LS220 EOS. In Bruenn et al. (2016) this model (as well as all others in that study) started shock expansion around 150 ms, while in Summa et al. (2016) this happened only around 300-350 ms.
Although the diversity of these results is quite considerable, this is not too astonishing in view of the fact that few-percent variations in a single neutrino interaction channel may already shift the time of shock runaway, , by several hundred milliseconds.
5 Summary and conclusions
In this study we used 1D (spherically symmetric) and 2D (axisymmetric) models to compare the relatively new Aenus-Alcar code (Just et al., 2015), which incorporates the fully multidimensional M1 approximation for neutrino transport, against the well-established Prometheus-Vertex code (Rampp & Janka, 2002; Buras et al., 2006), which employs an accurate Boltzmann solver restricted to the ray-by-ray+ (RbR+) approximation that neglects non-radial neutrino flux components. We compared with Vertex by mimicking the RbR+ approximation in Alcar and we tested the RbR+ approximation by comparing to the fully multidimensional version of Alcar. Moreover, we investigated the impact of other modeling variations in the neutrino transport that are frequently used by CCSN modelers, namely neglecting inelastic neutrino-electron scattering, simplifying pair-processes by assuming target neutrinos to be in isotropic equilibrium, applying strangeness and many-body corrections to neutrino-nucleon scattering, and ignoring velocity-dependent as well as gravitational redshift terms in the transport equations. Starting from spherically symmetric and non-rotating progenitor models, asymmetries are not expected to develop during the core-collapse phase, which therefore can be simulated by using computationally less demanding 1D calculations. Moreover, since the impact of some of our modeling variations on the (one-dimensional) core collapse has already been studied in detail before (e.g. Lentz et al., 2012b, a), we only focus here on the ramifications of these variations on the post-bounce evolution, by initializing all post-bounce models with the same collapse model. Finally, in order to obtain information about the degree of stochasticity in our comparison study we repeated some of our simulations several times starting with different initial random perturbation patterns (but the same perturbation amplitudes). With respect to the questions raised in the introduction, we obtained the following results:
In 1D the agreement between Alcar and Vertex is found to be excellent concerning nearly all features. The only noteworthy differences are a more energetic neutrino burst, slightly (%) higher luminosities, and a hotter and (by km) more compact proto-neutron star (PNS) in Alcar.
In the two examined 2D models of the SASI-dominated s20 progenitor with and without neutrino-electron scattering, the agreement found between Alcar in the RbR+ mode and Vertex remains very good concerning the neutrino emission, the PNS contraction, the heating conditions, and the explosion times. The two last mentioned features are subject to substantial stochastic scatter, the degree of which is, consistently in both codes, stronger for the exploding models with neutrino-electron scattering than for the non-exploding models without. Similar to the 1D case, the PNS radius in the Alcar models is again smaller throughout by km compared to the Vertex models. Although PNS convection has almost the same impact on the PNS radius and neutrino luminosities with Alcar and Vertex, the associated kinetic energies are higher by a factor of a few in Alcar than in Vertex. The origin of this discrepancy is not related to the RbR+ approximation and needs to be found in future work.
When comparing Alcar models with and without the RbR+ approximation, we could not observe any significant sensitivity of PNS convection and of the (angle-averaged) neutrino emission on the use of RbR+. Concerning the explosion behavior, we find that the RbR+ approximation only becomes noticeable once long-lived (polar) hot spots appear, above which material is heated more efficiently than without RbR+. In the investigated models hot spots are formed as a result of large-scale, low-order fluid modes that accumulate matter near the poles of the neutrinosphere. We observe clear explosion-promoting consequences for the SASI-dominated s20 models, but only a weak impact for the convection-dominated s9 models. We interpret this as a consequence of the tendency that linear sloshing modes may be more efficient than convection-driven modes in creating long-lasting accretion hot spots on the proto-neutron star surface. However, even in the SASI-dominated models the net impact of RbR+ remains manageable and quantitatively comparable to typical modeling variations in the microphysics sector (see next item).
Simplifying pair processes, including neutrino-electron scattering, as well as adopting the strangeness and many-body corrections during the post-bounce evolution all have, roughly in this order of relevance, a significant explosion-facilitating impact in our 2D models, while in 1D these modifications result in rather small changes of the shock radius. The net effect on the explosion times of combining the two considered opacity simplifications (regarding pair processes and neutrino-electron scattering) is smaller than that of using either simplification individually (cp. models s20-rbr, s20-rbr-nones, s20-rbr-pp, s20-rbr-pp-nones in Table 1). The timescale ratio, , in 1D predicts on a qualitative level remarkably well the impact of each modeling variation for the 2D models. In contrast, the shock radius can be shifted both to lower (“str” and “mb” models) or higher (“pp” models) values in 1D for microphysics variations that promote an explosion in 2D. This suggests that is a more powerful diagnostic quantity than the shock trajectory when estimating the impact of modeling variations using computationally less demanding 1D models.
Ignoring velocity-dependent terms and gravitational redshift in the transport equations during the post-bounce evolution reduces the neutrino fluxes noticeably and the rms-energies barely as measured in the comoving frame of infalling material in the gain region. In 2D these explosion-hampering features can, however, be (over-)compensated, because PNS contraction is accelerated owing to less vigorous PNS convection. The net effect on the explosion is case dependent and becomes more pessimistic for more massive PNSs, which imply higher velocities and stronger redshift: For an s20 model that originally explodes early when the PNS is still less massive, the explosion sets in earlier, while for an originally late exploding model the explosion lacks entirely. For the low-mass s9 models the aforementioned effects compensate and there is barely any visible impact on the shock trajectory.
The conditions in the gain layer (e.g. the lateral kinetic energy or timescale ratio) and the explosion times can be subject to substantial random variations of several hundred milliseconds, in agreement with the findings by Cardall & Budiardja (2015) based on simplified models. The amplitudes of temporal fluctuations and the scatter in the explosion times increase with the proximity to criticality for a given model, and they turned out to be much greater for the SASI-dominated s20 models than for the convection-dominated s9 models.
The result that our systematic comparison reveals good agreement between Alcar and Vertex is reassuring with respect to the numerical implementation of the codes and employed input physics, and it may also be considered as mutual support of the basic viability of the approximations made in both schemes, namely M1 and RbR+, at least for problems and setups similar to those considered in this study.
The tests of RbR+ confirm the previous suspicion (Dolence et al., 2015; Sumiyoshi et al., 2015; Skinner et al., 2016) that RbR+ can facilitate explosions, but they do so only partly, because our convection-dominated models are only marginally affected and a corresponding comparison in the more realistic three-dimensional case has yet to be conducted. Also, one should keep in mind that our comparison test of RbR+ bears uncertainties that are connected to the still incompletely known accuracy of fully multidimensional M1. For instance, it might be possible that M1 over- or underestimates the lateral dilution of radiation emitted from polar hot spots, in which case the observed differences between the SASI-dominated s20 models would presumably be over- or underrated here, respectively.
Remarkably, significant effects of RbR+ only seem to enter the evolution by means of long-lived, low-order multipole modes of the flow. In contrast, RbR+ effects triggered by temporary, stochastic downflows seem to remain weak on dynamical timescales (in agreement with the expectation of Buras et al., 2006), even though the instantaneous radiation field may exhibit significant local anisotropies with RbR+ (as found by Sumiyoshi et al., 2015). This suggests that in the more relevant 3D case the consequences of using RbR+ might be overall less dramatic than in 2D, because an artificial symmetry axis that fosters axis-parallel motions is absent. If confirmed, in turn, a reduced impact of RbR+ in 3D would speak in favor of existing 3D results obtained using the RbR+ approximation and would justify using RbR+ in the future.
Although the remaining results concerning the investigated simplifications are qualitatively consistent with previous studies emphasizing the importance of various aspects of neutrino transport (e.g. Buras et al., 2006; Müller et al., 2012a; Lentz et al., 2012b; Sumiyoshi et al., 2015; O’Connor & Couch, 2018; Burrows et al., 2018; Richers et al., 2017b; Kotake et al., 2018), the rather high sensitivity with respect to ostensibly small transport details combined with a considerable level of stochasticity is somewhat surprising. The following considerations are therefore not entirely new, but are certainly strengthened by this comparison study:
Profound knowledge of neutrino interaction rates is imperative. Further exploration of corrections to commonly used cross sections (e.g. Horowitz et al., 2017) and of new neutrino physics (e.g. Bollig et al., 2017) will more than likely have a significant leverage on future CCSN models. Our results suggest that a noticeable impact on the time of shock runaway can potentially be expected once a certain correction (or approximation) is large enough to cause a shift of the shock radius of just a few percent in corresponding spherically symmetric models.
The inclusion of (special and general) relativistic effects in multidimensional CCSN simulations is recommended, not only for obtaining a more precise time of shock runaway but also in order to achieve more reliable results for the PNS contraction and the luminosities and spectral properties of emitted neutrinos.
Careful analysis is necessary when comparing two models in the literature, particularly if these models use different approximations made in the neutrino sector and those differences are not rigorously tested and documented. In fact, such a comparison may turn out to be extremely difficult and quite amenable to premature conclusions. For instance, large differences in explosion times might be caused just by a different treatment of a certain neutrino interaction channel and may therefore not necessarily reflect more serious code differences. Likewise, good agreement between explosion times may not necessarily imply that all the individual differences are small; it could just happen that some differences accidentally cancel each other and the net impact is small. Moreover, due to their strongly non-linear nature, certain neutrino reactions may become more or less relevant when coupled with other reactions (see, e.g., Lentz et al., 2012b).
Tightly related to the previous item, the uncertainties due to stochasticity need to be estimated and accounted for whenever drawing conclusions regarding mutually different physics ingredients of two simulations in 2D. Obviously, performing expensive stochasticity tests is barely feasible (and might not be as important) for the currently most detailed 3D models. However, with the number of available models growing and the wealth of experience increasing also in 3D applications, more light will be shed on the quantitative impact of stochasticity and its dependence on the progenitor model and on other conditions.
Nonetheless, while our results certainly suggest to exercise caution when interpreting the runaway times of CCSN models using different neutrino treatments, we should also be aware that the level of sensitivity might not always be as dramatic as seen in the present study. The progenitor models chosen here, in particular the s20 model, might be more marginal than other, possibly more representative models. Moreover, once the set of physics ingredients leads to more robustly exploding models, details may tend to matter less, which means that errors introduced by some approximations (such as RbR+, M1, or neglecting neutrino-electron scattering) would become less significant.
A few final remarks are in order concerning the approximate M1 closure. Although Alcar with RbR+ and Vertex show very good agreement, at this point we are unable to quantify the error introduced in dynamical simulations by the fully multidimensional M1 closure (see, however, Richers et al., 2017b, for a comparison of M1 with different Boltzmann solvers based on stationary configuration snapshots). At least the error cannot be dramatically large in all cases of models and progenitors, since otherwise we would have seen stronger differences between the s9 models with and without RbR+. Nevertheless, it is known (e.g. Pons et al., 2000) that M1 can lead to unphysical features in the case of crossing radiation beams, and comparisons with more sophisticated transfer schemes in the case of black-hole torus systems (Appendix of Just et al., 2015a) or differentially rotating neutron stars (Foucart, 2018) as remnants of neutron-star mergers suggest that the accuracy of M1 systematically decreases with increasing geometric complexity of the source. The M1 scheme might therefore be somewhat more accurate for ordinary CCSNe than for more exotic ones with high rotation rates and/or strong magnetization, for which significant neutrino-emission asymmetries appear. Still, even for those cases the M1 method is currently one of the most attractive schemes regarding accuracy and computational efficiency, considering that for highly aspherical geometries the RbR+ approximation is not advisable and given the scarcity of alternative, computationally feasible multidimensional transport methods. In any case, comparisons with more accurate methods based on dynamical simulations are needed in order to identify the relevant shortcomings of M1, and to assess for which physics questions it will be necessary to employ more involved and expensive Boltzmann solvers.
More comparisons such as the present one are needed to understand differences in existing multidimensional CCSN results obtained worldwide. These comparisons should include more codes using different transport approximations (e.g. FLD, IDSA), as well as discretization schemes for the hydro (e.g. cartesian or cylindrical or spherical polar grid) and the transport (e.g. tangent-ray, discrete ordinate, Monte Carlo). Our results suggest that convection-dominated progenitors, which seem to come with an overall smaller level of stochasticity, might be more suited for such studies than SASI-dominated progenitors.
OJ is indebted to Thomas Ertl, Jérôme Guilet, Raphael Hix, Remi Kazeromi, Tobias Melson, Tomoya Takiwaki, Kohsuke Sumiyoshi, and Alexander Summa for helpful discussions. Moreover, we are grateful to Evan O’Connor and Almudena Arcones whose initiation of a community wide comparison project inspired parts of this work. The work of OJ, THJ, RB, and RG was supported by the European Research Council through grant ERC AdG 341157-COCO2CASA and the Cluster of Excellence âUniverseâ EXC 153, moreover for OJ and THJ by the Max-PlanckâPrinceton Center for Plasma Physics (MPPC). MO acknowledges support from the European Research Council (grant CAMAP-259276) and from the Spanish Ministry of Economy and Finance and the Valencian Community grants under grants AYA2015-66899-C2-1-P and PROMETEOII/2014-069. OJ and SN are supported by the Special Postdoctoral Researchers (SPDR) program and the iTHEMS cluster at RIKEN. We acknowledge computational support by the Max Planck Computing and Data Facility (MPCDF) and the HOKUSAI supercomputer at RIKEN.
- Arnett (1977) Arnett W. D., 1977, ApJ, 218, 815
- Audit et al. (2002) Audit E., Charrier P., Chièze J. ., Dubroca B., 2002, preprint (arXiv:astro-ph/0206281)
- Bethe & Wilson (1985) Bethe H. A., Wilson J. R., 1985, ApJ, 295, 14
- Blondin & Mezzacappa (2006) Blondin J. M., Mezzacappa A., 2006, ApJ, 642, 401
- Blondin et al. (2003) Blondin J. M., Mezzacappa A., DeMarino C., 2003, ApJ, 584, 971
- Bollig et al. (2017) Bollig R., Janka H.-T., Lohs A., Martínez-Pinedo G., Horowitz C. J., Melson T., 2017, Physical Review Letters, 119, 242702
- Bruenn (1985) Bruenn S. W., 1985, ApJS, 58, 771
- Bruenn & Mezzacappa (1997) Bruenn S. W., Mezzacappa A., 1997, Phys. Rev. D, 56, 7529
- Bruenn et al. (2016) Bruenn S. W., et al., 2016, ApJ, 818, 123
- Buras et al. (2006) Buras R., Rampp M., Janka H.-T., Kifonidis K., 2006, A&A, 447, 1049
- Burrows et al. (2000) Burrows A., Young T., Pinto P., Eastman R., Thompson T. A., 2000, ApJ, 539, 865
- Burrows et al. (2018) Burrows A., Vartanyan D., Dolence J. C., Skinner M. A., Radice D., 2018, Space Science Reviews, 214, 33
- Cabezón et al. (2018) Cabezón R. M., Pan K.-C., Liebendörfer M., Kuroda T., Ebinger K., Heinimann O., Thielemann F.-K., Perego A., 2018, preprint (arXiv:1806.09184)
- Cardall & Budiardja (2015) Cardall C. Y., Budiardja R. D., 2015, ApJL, 813, L6
- Cardall et al. (2013) Cardall C. Y., Endeve E., Mezzacappa A., 2013, Phys. Rev. D, 87, 103004
- Cernohorsky (1994) Cernohorsky J., 1994, ApJ, 433, 247
- Colella & Woodward (1984) Colella P., Woodward P. R., 1984, Journal of Computational Physics, 54, 174
- Colgate & White (1966) Colgate S. A., White R. H., 1966, ApJ, 143, 626
- Couch (2013) Couch S. M., 2013, ApJ, 775, 35
- Couch & Ott (2015) Couch S. M., Ott C. D., 2015, ApJ, 799, 5
- Couch et al. (2015) Couch S. M., Chatzopoulos E., Arnett W. D., Timmes F. X., 2015, ApJL, 808, L21
- Dessart et al. (2006) Dessart L., Burrows A., Livne E., Ott C. D., 2006, ApJ, 645, 534
- Dolence et al. (2015) Dolence J. C., Burrows A., Zhang W., 2015, ApJ, 800, 10
- Ertl et al. (2016) Ertl T., Janka H.-T., Woosley S. E., Sukhbold T., Ugliano M., 2016, ApJ, 818, 124
- Fernández (2015) Fernández R., 2015, MNRAS, 452, 2071
- Fernández et al. (2014) Fernández R., Müller B., Foglizzo T., Janka H.-T., 2014, MNRAS, 440, 2763
- Foglizzo (2002) Foglizzo T., 2002, A&A, 392, 353
- Foglizzo et al. (2006) Foglizzo T., Scheck L., Janka H.-T., 2006, ApJ, 652, 1436
- Foglizzo et al. (2007) Foglizzo T., Galletti P., Scheck L., Janka H.-T., 2007, ApJ, 654, 1006
- Foglizzo et al. (2015) Foglizzo T., et al., 2015, PASA, 32, e009
- Foucart (2018) Foucart F., 2018, MNRAS, 475, 4186
- Foucart et al. (2015) Foucart F., et al., 2015, Phys. Rev. D, 91, 124021
- Fryxell et al. (1989) Fryxell B., Müller E., W. A., 1989, Hydrodynamics and Nuclear Burning, MPA Green Report 449, Garching
- Guilet & Foglizzo (2012) Guilet J., Foglizzo T., 2012, MNRAS, 421, 546
- Hanke et al. (2012) Hanke F., Marek A., Müller B., Janka H.-T., 2012, ApJ, 755, 138
- Hannestad & Raffelt (1998) Hannestad S., Raffelt G., 1998, ApJ, 507, 339
- Horowitz (1997) Horowitz C. J., 1997, Phys. Rev. D, 55, 4577
- Horowitz (2002) Horowitz C. J., 2002, Phys. Rev. D, 65, 043001
- Horowitz et al. (2017) Horowitz C. J., Caballero O. L., Lin Z., O’Connor E., Schwenk A., 2017, Phys. Rev. C, 95, 025801
- Janka (2001) Janka H.-T., 2001, A&A, 368, 527
- Janka & Hillebrandt (1989) Janka H.-T., Hillebrandt W., 1989, A&AS, 78, 375
- Janka et al. (2016) Janka H.-T., Melson T., Summa A., 2016, Annual Review of Nuclear and Particle Science, 66, 341
- Just et al. (2015a) Just O., Bauswein A., Pulpillo R. A., Goriely S., Janka H.-T., 2015a, MNRAS, 448, 541
- Just et al. (2015) Just O., Obergaulinger M., Janka H.-T., 2015, MNRAS, 453, 3386
- Kazeroni et al. (2017) Kazeroni R., Guilet J., Foglizzo T., 2017, MNRAS, 471, 914
- Kotake et al. (2018) Kotake K., Takiwaki T., Fischer T., Nakamura K., Martínez-Pinedo G., 2018, The Astrophysical Journal, 853, 170
- Kraichnan (1967) Kraichnan R. H., 1967, Physics of Fluids, 10, 1417
- Kuroda et al. (2016) Kuroda T., Takiwaki T., Kotake K., 2016, ApJS, 222, 20
- Lattimer & Swesty (1991) Lattimer J. M., Swesty F., 1991, Nuclear Physics A, 535, 331
- Lentz et al. (2012a) Lentz E. J., Mezzacappa A., Messer O. E. B., Liebendörfer M., Hix W. R., Bruenn S. W., 2012a, ApJ, 747, 73
- Lentz et al. (2012b) Lentz E. J., Mezzacappa A., Messer O. E. B., Hix W. R., Bruenn S. W., 2012b, ApJ, 760, 94
- Lentz et al. (2015) Lentz E. J., et al., 2015, ApJL, 807, L31
- Levermore (1984) Levermore C. D., 1984, Journal of Quantitative Spectroscopy and Radiative Transfer, 31, 149
- Levermore & Pomraning (1981) Levermore C. D., Pomraning G. C., 1981, ApJ, 248, 321
- Liebendörfer et al. (2004) Liebendörfer M., Messer O. E. B., Mezzacappa A., Bruenn S. W., Cardall C. Y., Thielemann F.-K., 2004, ApJS, 150, 263
- Liebendörfer et al. (2005) Liebendörfer M., Rampp M., Janka H., Mezzacappa A., 2005, ApJ, 620, 840
- Liebendörfer et al. (2009) Liebendörfer M., Whitehouse S. C., Fischer T., 2009, ApJ, 698, 1174
- Marek et al. (2006) Marek A., Dimmelmeier H., Janka H.-T., Müller E., Buras R., 2006, A&A, 445, 273
- Marek et al. (2009) Marek A., Janka H.-T., Müller E., 2009, A&A, 496, 475
- Melson et al. (2015a) Melson T., Janka H.-T., Marek A., 2015a, ApJL, 801, L24
- Melson et al. (2015b) Melson T., Janka H.-T., Bollig R., Hanke F., Marek A., Müller B., 2015b, ApJL, 808, L42
- Mezzacappa & Bruenn (1993) Mezzacappa A., Bruenn S. W., 1993, ApJ, 405, 669
- Mignone (2014) Mignone A., 2014, Journal of Computational Physics, 270, 784
- Minerbo (1978) Minerbo G. N., 1978, J. Quant. Spec. Radiat. Transf., 20, 541
- Müller (2016) Müller B., 2016, PASA, 33, e048
- Müller & Janka (2015) Müller B., Janka H.-T., 2015, MNRAS, 448, 2141
- Müller et al. (2012a) Müller B., Janka H.-T., Marek A., 2012a, ApJ, 756, 84
- Müller et al. (2012b) Müller B., Janka H.-T., Heger A., 2012b, ApJ, 761, 72
- Müller et al. (2016) Müller B., Viallet M., Heger A., Janka H.-T., 2016, ApJ, 833, 124
- Müller et al. (2017) Müller B., Melson T., Heger A., Janka H.-T., 2017, MNRAS, 472, 491
- Nagakura et al. (2018) Nagakura H., et al., 2018, ApJ, 854, 136
- Nakamura et al. (2015) Nakamura K., Takiwaki T., Kuroda T., Kotake K., 2015, PASJ, 67, 107
- O’Connor (2015) O’Connor E., 2015, ApJS, 219, 24
- O’Connor & Couch (2018) O’Connor E. P., Couch S. M., 2018, ApJ, 854, 63
- O’Connor & Ott (2010) O’Connor E., Ott C. D., 2010, Classical and Quantum Gravity, 27, 114103
- O’Connor et al. (2017) O’Connor E., Horowitz C. J., Lin Z., Couch S., 2017, in Marcowith A., Renaud M., Dubner G., Ray A., Bykov A., eds, IAU Symposium Vol. 331, Supernova 1987A:30 years later - Cosmic Rays and Nuclei from Supernovae and their Aftermaths. pp 107–112, arXiv:1712.08253, doi:10.1017/S1743921317004586
- O’Connor et al. (2018) O’Connor E., et al., 2018, preprint (arXiv:1806.041750)
- Obergaulinger (2008) Obergaulinger M., 2008, Dissertation, Technische Universität München, München
- Obergaulinger & Aloy (2017) Obergaulinger M., Aloy M. Á., 2017, MNRAS, 469, L43
- Ono et al. (2013) Ono M., Nagataki S., Ito H., Lee S.-H., Mao J., Hashimoto M.-a., Tolstov A., 2013, ApJ, 773, 161
- Ott et al. (2017) Ott C. D., Roberts L. F., da Silva Schneider A., Fedrow J. M., Haas R., Schnetter E., 2017, preprint (arXiv:1712.01304)
- Pan et al. (2018) Pan K.-C., Mattes C., O’Connor E. P., Couch S. M., Perego A., Arcones A., 2018, preprint (arXiv:1806.10030)
- Perego et al. (2016) Perego A., Cabezón R. M., Käppeli R., 2016, ApJS, 223, 22
- Pons et al. (1998) Pons J. A., Miralles J. A., Ibanez J. M. A., 1998, A&AS, 129, 343
- Pons et al. (2000) Pons J. A., Ibáñez J. M., Miralles J. A., 2000, MNRAS, 317, 550
- Radice et al. (2018) Radice D., Abdikamalov E., Ott C. D., Mösta P., Couch S. M., Roberts L. F., 2018, Journal of Physics G Nuclear Physics, 45, 053003
- Radice et al. (2017) Radice D., Abdikamalov E., Ott C. D., Moesta P., Couch S. M., Roberts L. F., 2017, preprint (arXiv:1710.01282)
- Rampp & Janka (2002) Rampp M., Janka H., 2002, A&A, 396, 361
- Richers et al. (2017a) Richers S., Ott C. D., Abdikamalov E., O’Connor E., Sullivan C., 2017a, Phys. Rev. D, 95, 063019
- Richers et al. (2017b) Richers S., Nagakura H., Ott C. D., Dolence J., Sumiyoshi K., Yamada S., 2017b, ApJ, 847, 133
- Roberts et al. (2016) Roberts L. F., Ott C. D., Haas R., O’Connor E. P., Diener P., Schnetter E., 2016, ApJ, 831, 98
- Ruffert et al. (1996) Ruffert M., Janka H.-T., Schaefer G., 1996, A&A, 311, 532
- Scheck et al. (2008) Scheck L., Janka H.-T., Foglizzo T., Kifonidis K., 2008, A&A, 477, 931
- Sekiguchi et al. (2015) Sekiguchi Y., Kiuchi K., Kyutoku K., Shibata M., 2015, Phys. Rev. D, 91, 064059
- Shibata et al. (2011) Shibata M., Suwa Y., Kiuchi K., Ioka K., 2011, ApJL, 734, L36
- Skinner et al. (2016) Skinner M. A., Burrows A., Dolence J. C., 2016, ApJ, 831, 81
- Steiner et al. (2013) Steiner A. W., Hempel M., Fischer T., 2013, ApJ, 774, 17
- Sukhbold et al. (2016) Sukhbold T., Ertl T., Woosley S. E., Brown J. M., Janka H.-T., 2016, ApJ, 821, 38
- Sumiyoshi et al. (2015) Sumiyoshi K., Takiwaki T., Matsufuru H., Yamada S., 2015, ApJS, 216, 5
- Summa et al. (2016) Summa A., Hanke F., Janka H.-T., Melson T., Marek A., Müller B., 2016, ApJ, 825, 6
- Suwa et al. (2016) Suwa Y., Yamada S., Takiwaki T., Kotake K., 2016, ApJ, 816, 43
- Takiwaki et al. (2014) Takiwaki T., Kotake K., Suwa Y., 2014, ApJ, 786, 83
- Tamborra et al. (2017) Tamborra I., Hüdepohl L., Raffelt G. G., Janka H.-T., 2017, ApJ, 839, 132
- Thompson et al. (2003) Thompson T. A., Burrows A., Pinto P. A., 2003, ApJ, 592, 434
- Toro (1997) Toro E., 1997, Riemann solvers and numerical methods for fluid dynamics: a practical introduction. Springer Verlag
- Wongwathanarat et al. (2017) Wongwathanarat A., Janka H.-T., Müller E., Pllumbi E., Wanajo S., 2017, ApJ, 842, 13
- Woosley & Heger (2007) Woosley S. E., Heger A., 2007, Physics Reports, 442, 269
- Woosley & Heger (2015) Woosley S. E., Heger A., 2015, ApJ, 810, 34
- Yamada et al. (1999) Yamada S., Janka H.-T., Suzuki H., 1999, A&A, 344, 533
- Yueh & Buchler (1977) Yueh W. R., Buchler J. R., 1977, ApJ, 217, 565
Appendix A Treatment of source terms in Alcar
In Alcar the overall time integration of the neutrino moments is performed together with the hydrodynamic variables using a 2nd-order Runge-Kutta scheme, in which all terms except some neutrino-interaction rates are treated explicitly in time. For the motivation of the overall integration method we refer the reader to Just et al. (2015). Here we describe the specific treatment of neutrino-interaction rates employed for the simulations in this paper.
Ignoring for now the Runge-Kutta sub-stepping, we advance the hydrodynamic equations, Eqs. (1), and neutrino-transport equations, Eqs. (3), from an old time step, , to a new (partial) time step, , as follows:
where (with and denoting the neutrino species and energy group, respectively) is the vector of hydrodynamic and neutrino-transport variables at the old/new time step, are the neutrino-source terms, which may depend on both old and new variables, and represents all remaining terms, which depend only on the old variables. Now if in Eq. (16) all neutrino-interaction terms for the three species, four moments, and energy groups were to be treated fully implicitly, at each spatial grid point a non-linear system of equations of rank (the additional two variables are and ) would need to be solved, which involves the inversion of a non-sparse Jacobian possibly multiple times per integration step. Given that the time steps used in Alcar are rather short (because of the explicit integration of the non-local divergence terms) and the total number of integration steps is therefore large, the computational cost would in this case soon become very large. This is why we avoid an implicit integration of energy-bin coupling source terms wherever justified, in the manner that is described below.
For all emission/absorption as well as iso-energetic scattering processes the source terms entering the moment equations, Eqs. (3), are given by:
where is the equilibrium energy density corresponding to a Fermi-Dirac distribution, is the sum of all absorption opacities corrected for stimulated absorption, and is the sum of opacities for iso-energetic scattering. We follow Rampp & Janka (2002) for the computation of these opacities, except that we additionally account for weak magnetism and nucleon recoil (Horowitz, 2002) and, in selected models (cf. Sec. 2.3), for strangeness corrections and many-body corrections in a way described in (Horowitz et al., 2017). For neutral-current reactions of neutrinos (representing the four heavy-lepton neutrinos) we use as effective opacity the arithmetic average of the opacity of a heavy-lepton neutrino and that of its antiparticle. In computing the source terms, Eqs. (17), the neutrino energy- and flux-densities are treated implicitly and the opacities explicitly. The equilibrium energy density, , is usually treated explicitly, i.e. using old values of and . Occasionally, however, in regions with strong neutrino-matter coupling and short fluid-dynamical timescales the difference between equilibrium energies at the old and new time step becomes so large that numerical oscillations would result for an explicit treatment of . In these (rare) cases, which we detect using a criterion based on the density and the relative change of gas energy due to neutrino interactions, we perform an additional intermediate step to obtain improved values of and (and therefore ). In this step we solve Eqs. (16) for and using implicit equilibrium energies and under the simplifying assumption that the neutrino fluxes vanish.
Inelastic scattering of (all types of) neutrinos off electrons and positrons is implemented as in Yueh & Buchler (1977); Bruenn (1985); Rampp & Janka (2002), i.e. using a Legendre expansion in the scattering angle of neutrinos up to 1st order. For the multidimensional generalization of the formalism of Yueh & Buchler (1977); Bruenn (1985); Rampp & Janka (2002) (who assumed spherical symmetry and therefore vanishing non-radial flux-vector components and non-diagonal Eddington-tensor components) we use the expressions given in Cernohorsky (1994). We obtain the Legendre coefficients of up to 1st order for each initial-state and final-state neutrino energy and each neutrino species by table interpolation. The scattering rates are computed explicitly, i.e. using Legendre coefficients and neutrino moments from the old time step. In order for the explicit integration to remain stable, the rates for inelastic neutrino-lepton scattering are reduced at high densities (where these rates are subdominant compared to charged-current reactions for conditions considered in this paper) by a factor of (following the suggestion by O’Connor, 2015).
Finally, for pair processes we again follow Rampp & Janka (2002) as closely as possible121212We note the following misprints in the Appendix of Rampp & Janka (2002), which are correctly accounted for in both simulation codes, Alcar and Vertex: In Eq.(A.20), i.e. the multipole expansion of the 1st-moment source term for pair-processes, the factor must be replaced by 1, in Eq.(A.39) the quantity erroneously denoted as production coefficient, , is actually the annihilation coefficient, , and in Eq.(A.41) the factor must be replaced by 1., where electron-positron annihilation is implemented based on Pons et al. (1998) and nucleon-nucleon bremsstrahlung is implemented based on Hannestad & Raffelt (1998). We expand the interaction kernel in , where n and are the propagation unit vectors of the considered neutrino and its annihilation partner, respectively. The resulting expansion of the pair-process source terms for the moments and of a given neutrino species at energy reads (up to 2nd order in ):