Gray two-moment neutrino transport: Comprehensive tests and improvements for supernova simulations

Haakon Andresen; Evan P. O’Connor; Oliver Eggenberger Andersen; Sean M. Couch

doi:10.1051/0004-6361/202449776

Home

All issues

Volume 687 (July 2024)

A&A, 687 (2024) A55

Full HTML

Open Access

Issue		A&A Volume 687, July 2024


Article Number		A55
Number of page(s)		18
Section		Numerical methods and codes
DOI		https://doi.org/10.1051/0004-6361/202449776
Published online		28 June 2024

A&A, 687, A55 (2024)

Gray two-moment neutrino transport: Comprehensive tests and improvements for supernova simulations

Haakon Andresen¹, Evan P. O’Connor¹, Oliver Eggenberger Andersen¹ and Sean M. Couch²^,3^,4

¹ The Oskar Klein Centre, Department of Astronomy Stockholm University, AlbaNova, 106 91 Stockholm, Sweden
e-mail: haakon.andresen@astro.su.se
² Department of Physics and Astronomy, Michigan State University, East Lansing, MI 48824, USA
³ Department of Computational Mathematics, Science, and Engineering, Michigan State University, East Lansing, MI 48824, USA
⁴ Facility for Rare Isotope Beams, Michigan State University, East Lansing, MI 48824, USA

Received: 28 February 2024
Accepted: 10 April 2024

Abstract

Aims. In this work we extended an energy-integrated neutrino transport method to facilitate efficient, yet precise, modeling of compact astrophysical objects. We particularly focus on core-collapse supernovae.

Methods. We implemented a gray neutrino-transport framework from the literature into FLASH and performed a detailed evaluation of its accuracy in core-collapse supernova simulations. Based on comparisons with results from simulations using energy-dependent neutrino transport, we incorporated several improvements to the original scheme.

Results. Our analysis shows that our gray neutrino transport method successfully reproduces key aspects from more complex energy-dependent transport across a variety of progenitors and equations of state. We find both qualitative and reasonable quantitative agreement with multi-group M1 transport simulations. However, the gray scheme tends to slightly favor shock revival. In terms of gravitational wave and neutrino signals, there is a good alignment with the energy-dependent transport, although we find 15–30% discrepancies in the average energy and luminosity of heavy-lepton neutrinos. Simulations using the gray transport are around four times faster than those using energy-dependent transport.

Key words: gravitational waves / hydrodynamics / neutrinos / radiative transfer / supernovae: general

© The Authors 2024

Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 Introduction

Core-collapse supernovae are important for several facets of astrophysics: they play a key role in the chemical (Edmunds 2017; Thielemann et al. 2018) and dynamical evolution (Smith et al. 2018; Bacchini et al. 2023) of galaxies and they are the progenitors of every binary system observed to date by the LIGO-Virgo-KAGRA collaboration (Abbott et al. 2016b,c,a, 2017d,a,b,c, 2019, 2020b,c,a,d, 2021b,a,a, 2023, 2024). To thoroughly understand core-collapse supernovae and their impact in astrophysics, it is essential to systematically examine the varied outcomes that occur in the final stages of the evolution of massive stars.

In the last decade, the efforts of the modeling community have led to several successful supernova explosions in three-dimensional simulations with state-of-the-art microphysics (Melson et al. 2015a,b; Roberts et al. 2016; Summa et al. 2018; Vartanyan et al. 2019, 2022; Burrows et al. 2019, 2020). Despite recent advances, there are still several open questions and challenges left to tackle in the core-collapse field. The ultimate goal is to simulate a star from the onset of iron-core collapse, through shock revival, until the shock breaks out of the star and self-consistently predicting every aspect of the observables. Reaching this goal will require an improved understanding of the stellar progenitors (Arnett & Meakin 2011; Couch & Ott 2015; Cristini et al. 2017; Jones et al. 2017; Müller et al. 2016), understanding the importance of neutrino oscillations (Izaguirre et al. 2017; Chakraborty et al. 2016; Capozzi et al. 2019; Johns et al. 2020; Chakraborty & Chakraborty 2020; Bhattacharyya & Dasgupta 2021; Capozzi et al. 2020, 2022; Martin et al. 2021; Johns 2023; Xiong et al. 2023a,b; Richers & Sen 2022; Dedin Neto et al. 2023; Liu et al. 2023; Cornelius et al. 2024; Shalgar & Tamborra 2024; Ehring et al. 2023; Akaho et al. 2024), accurate microphysics (Sumiyoshi et al. 2005; Hempel et al. 2012; Fischer et al. 2014; Oertel et al. 2017; da Silva Schneider et al. 2020; Pascal et al. 2022; Suleiman et al. 2023), understanding the impact of rotation and magnetic fields (Kotake et al. 2011; Takiwaki & Kotake 2018; Summa et al. 2018; Jardine et al. 2022; Obergaulinger & Aloy 2021; Bugli et al. 2021; Reichert et al. 2023; Bugli et al. 2023; Buellet et al. 2023), and the ability to accurately simulate all the above-mentioned aspects. The input physics underlying core-collapse supernovae are complex and are subject to active research. Consequently, the underpinning of supernova simulations comes with inherent uncertainty.

On the one hand, state-of-the-art simulations are necessary to move the field forward, and the inclusion of ever more detailed physics in numerical simulations is sure to yield interesting results in the future. On the other hand, the most complex simulations are ill-suited for systematically studying the effects of the uncertainties in the input physics due to the high computational cost. First among the computationally expensive parts of core-collapse simulations is the transport of neutrinos from the optically thick regions of the forming neutron star to the outer layers of the core, from where they freestream out of the progenitor. Simplified neutrino transport methods have the potential to reduce the computational cost of numerical simulations significantly, but accurately tracking the radiation fields is important since neutrinos play a crucial role in the dynamics of core-collapse supernovae. Significant work has therefore been invested into improving the neutrino treatment in supernova simulations. Currently, the state-of-the-art consists of approximate solutions of the energy-dependent Boltzmann equation. Common approaches include moment schemes (O’Connor & Couch 2018b; Melson et al. 2015a; Just et al. 2015; Summa et al. 2018; Vartanyan et al. 2018; Kuroda et al. 2018, 2022; Just et al. 2018; Burrows et al. 2020; Wang & Burrows 2024), flux limited diffusion (Bruenn et al. 2020; Rahman et al. 2019, 2022), fast multi-group transport (Müller & Janka 2015), an improved leakage-equilibration-absorption scheme (Ardevol-Pulpillo et al. 2019), and the isotropic diffusion source approximation (Takiwaki et al. 2014; Pan et al. 2016; Kotake et al. 2018; Matsumoto et al. 2022). Monte Carlo methods (Abdikamalov et al. 2012; Richers et al. 2017; Kato et al. 2020) and solving the full Boltzmann equation (Liebendörfer et al. 2004; Nagakura et al. 2014, 2018, 2019; Iwakami et al. 2020; Harada et al. 2020; Akaho et al. 2021, 2023) has been attempted. However, due to the numerical cost and complexity, neither full Boltz-mann methods nor Monte Carlo methods have been widely adopted in simulations (see Mezzacappa et al. 2020 for a detailed review). Energy-integrated, or gray, schemes are often employed in simulations of binary neutron star mergers (Foucart et al. 2015, 2016; Fujibayashi et al. 2020, 2023; Radice et al. 2022; Hayashi et al. 2022; Zappa et al. 2023; Kiuchi et al. 2023, 2024; Curtis et al. 2024; Radice & Bernuzzi 2023; Schianchi et al. 2024; Foucart 2023). While gray transport methods were used for core-collapse supernovae in the past (Burrows et al. 2000; Fryer & Warren 2004; Scheck et al. 2006), they have largely been abandoned in the continuous effort to improve the neutrino treatment since it was seen as key to solving the supernova problem (producing successful explosions in numerical simulations).

Gray transport is significantly less computationally expensive than energy-dependent neutrino transport, which is an advantage if one wants to perform systematic studies. In this work we develop a gray moment scheme, based on methods commonly used in the binary merger literature (Foucart et al. 2015, 2016), and compare our results to fully energy-dependent neutrino transport. By comparing our gray scheme to a more sophisticated method, we were able to improve several aspects of the original gray transport and achieve good overall agreement between the two approaches. We started by implementing the energy-integrated neutrino transport suggested by Foucart et al. (2015, 2016) into FLASH. We then tested its accuracy, and suggest specific changes based on our test results.

This paper is organized as follows. We first briefly describe the current version of FLASH in Sect. 2. In Sect. 3, we summarize the neutrino transport problem and the method described by Foucart et al. (2016). We present numerical tests of our implementation of the neutrino transport scheme from Foucart et al. (2016) in Sect. 4. In Sect. 5, we present the changes we propose to the original transport scheme. Toward the end of Sect. 5 we present the same tests we initially performed for the scheme Foucart et al. (2016) for our modified method. We then describe the results from several two-dimensional (2D) test simulations in Sect. 6. In Sect. 7, we present the gravitational waves (GWs) and neutrino signals from a simulation utilizing our gray scheme and compare the results to signals from a corresponding simulation with energy-dependent transport. In Sect. 8, we discuss the performance of the gray transport, compared to the energy-dependent transport. Finally, in Sect. 9 we summarize our findings and present our conclusions.

2 FLASH

Our work builds upon the energy-dependent two-moment M1 neutrino transport (Cardall et al. 2013; Shibata et al. 2011) implemented in the FLASH code (O’Connor & Couch 2018b,a). We utilized a modified version of the FLASH framework (version 4; Fryxell et al. 2000), was specifically adapted for core-collapse simulations (Couch 2013; Couch & O’Connor 2014; O’Connor & Couch 2018a). FLASH implements adaptive mesh refinement, solves the Newtonian hydrodynamic equations using a modified general relativistic effective potential (case A in Marek et al. 2006), and evolves three neutrino species: electron neutrinos, electron anti-neutrinos, and a third species collecting all the heavy-lepton neutrinos into one species. Neutrino opacities were generated with the neutrino opacity library NuLib (O’Connor 2015; Sullivan et al. 2016). Heavy-lepton neutrinos are mainly produced in the proto-neutron star (PNS) through pair-production processes, whose rates depend on the distribution function of both neutrinos and anti-neutrinos. Since FLASH evolves a single species for the heavy-lepton neutrinos, the neutrino pair-process reactions are calculated in an approximate way following Burrows et al. (2006). An effective neutrino emission is calculated assuming isotropic emission, no final state blocking, and integrating over the energy of the anti-neutrino. The effective absorption is calculated through Kirchhoff’s law from the effective emissivity. We refer to O’Connor (2015) for details regarding the opacity calculations.

3 Neutrino transport

In the limit of massless neutrinos, the neutrino transport is described by the Boltzmann equation (Lindquist 1966), $p^{α} [\frac{\partial f_{v}}{\partial x^{α}} - Γ_{α γ}^{β} p^{γ} \frac{\partial f_{v}}{\partial p^{β}}] = {[\frac{\partial f_{v}}{\partial τ}]}_{coll} .$ ${p^\alpha }\left[ {{{\partial {f_v}} \over {\partial {x^\alpha }}} - \Gamma _{\alpha \gamma }^\beta {p^\gamma }{{\partial {f_v}} \over {\partial {p^\beta }}}} \right] = {\left[ {{{\partial {f_v}} \over {\partial \tau }}} \right]_{{\rm{coll}}}}.$ (1)

Here f = f(x^α, p^α) is the distribution function of the neutrinos, with x^α = (t, xⁱ) and p^α representing the four-position and the four-momentum of the neutrinos, respectively. $Γ_{α γ}^{β}$ $\Gamma _{\alpha \gamma }^\beta$ are the Christoffel symbols. The term on the right-hand side of Eq. (1) represents the scattering, absorption and emission of neutrinos (collectively referred to as collision processes). Note the subscript v represents the different neutrino species, in this work we evolve three neutrino species: electron neutrinos (v_e), electron anti-neutrinos $({\bar{v}}_{e})$ $\left( {{{\bar v}_e}} \right)$ , and a third species (v_x) representing the heavy-lepton neutrinos.

The complete transport problem spans a seven-dimensional phase space, which makes the problem prohibitively expensive to solve numerically. To circumvent this issue, we employ the truncated-moment formalism approximation (Thorne 1980). This formalism simplifies the problem by only evolving the lowest order moments of the Boltzmann equation. The problem is further simplified by integrating over energy and evolving the energy-integrated, or gray, moments.

3.1 Evolution equations

We summarize the method proposed by Foucart et al. (2016), but we refer to the original work for additional details and to O’Connor & Couch (2018a) for details on the M1 neutrino transport implemented in FLASH. Unlike Foucart et al. (2016), we operate within the Newtonian framework. Therefore, our formulation differs slightly from that of Foucart et al. (2016). In Sect. 5, we describe the changes we made to the original scheme by Foucart et al. (2016), but for clarity and to provide a summary of our starting point, we describe the original method below.

By evolving the number density of the neutrinos, in addition to the energy and momentum densities, it is possible to get a local estimate of the average energy of the neutrinos, which facilitates an approximation of the spectral shape of the neutrino spectrum (Foucart et al. 2016).

The evolution equations for the energy density (E), momentum density (F_i), and number density (N) are in the laboratory frame, $\partial_{t} E + \partial_{i} [α F^{i}] + α F^{i} \partial_{i} Φ = α (W [η - κ_{a} J] - κ_{H} H^{t}),$ ${\partial _t}E + {\partial _i}\left[ {\alpha {F^i}} \right] + \alpha {F^i}{\partial _i}\Phi = \alpha \left( {W\left[ {\eta - {\kappa _a}J} \right] - {\kappa _H}{H^t}} \right),$ (2) $\partial_{t} F^{i} + \partial_{j} [α P^{i j}] + α E \partial_{i} Φ = α (W [η - κ_{a} J] v^{i} - κ_{H} H^{i}),$ ${\partial _t}{F^i} + {\partial _j}\left[ {\alpha {P^{ij}}} \right] + \alpha E{\partial _i}\Phi = \alpha \left( {W\left[ {\eta - {\kappa _a}J} \right]{v^i} - {\kappa _H}{H^i}} \right),$ (3) $\partial_{t} N + \partial_{j} [α F_{N}^{j}] = α (η_{N} - \frac{κ_{N} J N}{W (E - F_{i} v^{i})}) .$ ${\partial _t}N + {\partial _j}\left[ {\alpha F_N^j} \right] = \alpha \left( {{\eta _N} - {{{\kappa _N}JN} \over {W\left( {E - {F_i}{v^i}} \right)}}} \right).$ (4)

Here η represents the energy-integrated emissivity, η_N is the energy-integrated number emissivity, and κ_N is the energy-averaged number absorption. κ_H denotes the sum of the energy-averaged absorption (κ_a) and energy-averaged scattering opacities (κ_s), κ_H = κ_a + κ_s. The energy-averaged opacities are detailed in Sect. 3.3. The fluid velocity is represented by υⁱ, W is the Lorentz factor, the gravitational potential is denoted with Φ, and the lapse is α = exp(Φ). The terms J and H^µ denote the energy density and momentum density in the fluid rest frame, respectively, and are related to the laboratory variables as follows: $J = W^{2} [E - 2 F^{i} v_{i} + v^{i} v^{j} P_{i j}],$ $J = {W^2}\left[ {E - 2{F^i}{v_i} + {v^i}{v^j}{P_{ij}}} \right],$ (5) $H^{t} = W^{3} [- (E - F^{i} v_{i}) v^{2} + F^{i} v_{i} - v_{i} v_{j} P^{i j}],$ ${H^t} = {W^3}\left[ { - \left( {E - {F^i}{v_i}} \right){v^2} + {F^i}{v_i} - {v_i}{v_j}{P^{ij}}} \right],$ (6) $H^{i} = W^{3} [- (E - F^{j} v_{j}) v^{i} + F^{i} - v_{j} P^{i j}] .$ ${H^i} = {W^3}\left[ { - \left( {E - {F^j}{v_j}} \right){v^i} + {F^i} - {v_j}{P^{ij}}} \right].$ (7)

In the above equations, P^ij are the second-order moments of the multi-pole expansion, which in the two-moment approach is determined by an analytic closure relation (O’Connor & Couch 2018a).

In Eq. (4), $F_{N}^{j}$ $F_N^j$ is the number momentum density that we assume takes the form of $F_{N}^{j} = \frac{J W v^{j}}{〈 ε 〉} + \frac{H^{j}}{〈 ε^{F} 〉},$ $F_N^j = {{JW{v^j}} \over {\langle \varepsilon \rangle }} + {{{H^j}} \over {\left\langle {{\varepsilon ^F}} \right\rangle }},$ (8)

where 〈ε〉 is the average neutrino energy and 〈ε^F〉 is the average energy of the neutrino flux (see Sect. 3.4).

3.2 Neutrino spectrum

The gray scheme of Foucart et al. (2016) is underpinned by the assumption that the neutrinos follow a Fermi-Dirac distribution $f_{v} (ε) = \frac{1}{1 + \exp ((ε - μ_{v}) / T_{v})},$ ${f_v}(\varepsilon ) = {1 \over {1 + \exp \left( {\left( {\varepsilon - {\mu _v}} \right)/{T_v}} \right)}},$ (9)

where ε represents the neutrino energy, µ_ν is the chemical potential of the neutrinos and is the T_ν the neutrino temperature. With the corresponding blackbody intensity $B_{v} (ε) \propto \frac{ε^{3}}{1 + \exp ((ε - μ_{v}) / T_{v})} .$ ${B_v}(\varepsilon ) \propto {{{\varepsilon ^3}} \over {1 + \exp \left( {\left( {\varepsilon - {\mu _v}} \right)/{T_v}} \right)}}.$ (10)

Accurate knowledge of the neutrino temperature is necessary for accurately describing the neutrino spectrum. One approach is to assume that the neutrino temperature is equal to the fluid temperature, but this can lead to errors in regimes where the neutrinos are not in equilibrium with the fluid.

By evolving the number density it is possible to get a better estimate of the neutrino temperature through the average energy of the neutrinos (Foucart et al. 2016), given by $〈 ε 〉 = W \frac{E - F_{i} v^{i}}{N} .$ $\langle \varepsilon \rangle = W{{E - {F_i}{v^i}} \over N}.$ (11)

For a blackbody spectrum, the average neutrino energy can be expressed in terms of the neutrino temperature and the ratio of the third-, and second-order Fermi integrals $〈 ε 〉 = \frac{F_{3} (η_{v})}{F_{2} (η_{v})} T_{v},$ $\langle \varepsilon \rangle = {{{F_3}\left( {{\eta _v}} \right)} \over {{F_2}\left( {{\eta _v}} \right)}}{T_v},$ (12)

where η_ν = µ_ν/T (not to be confused with the emissivity or number emissivity, η or η_N). The Fermi integrals are defined as follows¹: $F_{k} (η_{v}) = \int_{0}^{\infty} \frac{x^{k}}{1 + e^{x - η_{v}}} d x .$ ${F_k}\left( {{\eta _v}} \right) = \int_0^\infty {{{{x^k}} \over {1 + {{\rm{e}}^{x - {\eta _v}}}}}} {\rm{d}}x.$ (13)

We can then find the neutrino temperature by combining equations Eqs. (11) and (12). In our implementation of the gray scheme, we follow Foucart et al. (2015) and chose the following expression for the chemical potential $μ_{v} = μ_{v}^{eq} (1 - e^{- τ}),$ ${\mu _v} = \mu _v^{{\rm{eq}}}\left( {1 - {{\rm{e}}^{ - \tau }}} \right),$ (14)

where $μ_{γ}^{eq}$ $\mu _\gamma ^{{\rm{eq}}}$ is the chemical potential of the neutrinos assuming equilibrium with the fluid and τ is the optical depth of the neutrinos.

3.3 Neutrino opacities

The energy-integrated neutrino-matter source terms in Eqs. (2)– (4) play an important role in numerical simulations. However, the best choice for these quantities is not immediately obvious. The energy-integrated emissivity, which represents the total amount of neutrinos emitted by the fluid, can be defined as $η = \int η (ρ, Y_{e}, T, ε) d ε .$ $\eta = \int \eta \left( {\rho ,{Y_e},T,\varepsilon } \right){\rm{d}}\varepsilon .$ (15)

In the above equation, ρ, Y_e, and T represents the density, electron fraction, and temperature of the fluid, respectively. The energy-averaged opacities provide a more complex problem, because the energy average should be weighed by the neutrino spectrum. In other words, in Eqs. (2), and (3), $κ_{a} = \frac{\int κ_{a} (ε) J (ε) d ε}{\int J (ε) d ε},$ ${\kappa _a} = {{\int {{\kappa _a}} (\varepsilon )J(\varepsilon ){\rm{d}}\varepsilon } \over {\int J (\varepsilon ){\rm{d}}\varepsilon }},$ (16)

and $κ_{H} = \frac{\int κ_{H} (ε) H^{r} (ε) d ε}{\int H^{r} (ε) d ε} .$ ${\kappa _H} = {{\int {{\kappa _H}} (\varepsilon ){H^r}(\varepsilon ){\rm{d}}\varepsilon } \over {\int {{H^r}} (\varepsilon ){\rm{d}}\varepsilon }}.$ (17)

Here H^r is the radial component of the momentum density, which we use because our comparison of the energy-averaged opacities is based on 1D simulations. Consequently, even under the blackbody assumption, it is necessary to know the properties of the radiation field in order to perform the averaging. One approach is to assume that the neutrinos are in equilibrium with the fluid, which means that J ∝ B_ν, and define $κ_{a / H}^{eq} = \frac{\int B_{v} (T, μ_{v}, ε) κ_{a / H} (ρ, Y_{e}, T, ε) d ε}{\int B_{v} (T, μ_{v}, ε) d ε} .$ $\kappa _{a/H}^{{\rm{eq}}} = {{\int {{B_v}} \left( {T,{\mu _v},\varepsilon } \right){\kappa _{a/H}}\left( {\rho ,{Y_e},T,\varepsilon } \right){\rm{d}}\varepsilon } \over {\int {{B_v}} \left( {T,{\mu _v},\varepsilon } \right){\rm{d}}\varepsilon }}.$ (18)

The equilibrium assumption is denoted by the superscript “eq”. Assuming that the neutrinos are in equilibrium with the fluid works well in the optically thick regions, but less so in the semi-transparent and optically thin regions where the neutrino temperature can significantly differ from the fluid temperature. To account for the difference between the neutrino and fluid temperature, Foucart et al. (2016) proposed to weigh the opacities with the square of the ratio of the two temperatures and set $κ_{a / H} = κ_{a / H}^{eq} \frac{T_{v}^{2}}{T^{2}} .$ ${\kappa _{a/H}} = \kappa _{a/H}^{{\rm{eq}}}{{T_v^2} \over {{T^2}}}.$ (19)

The correction factor is based on the fact that the opacities typically go as energy squared, in other words, κ ~ ε². Defining the opacities in terms of Eq. (19) ensures that $κ_{a} = η J$ ${\kappa _a} = \eta J$ (20)

in the optically thick regions where the fluid is in equilibrium with the fluid (T_ν = T) and consequently guarantees detailed balance.

For the number emissivity and number absorption Foucart et al. (2016) chose $κ_{N} = κ_{a} \frac{η_{N}}{η} \frac{F_{3} (η_{v}) T}{F_{2} (η_{v})},$ ${\kappa _N} = {\kappa _a}{{{\eta _N}} \over \eta }{{{F_3}\left( {{\eta _v}} \right)T} \over {{F_2}\left( {{\eta _v}} \right)}},$ (21)

with κ_a from Eq. (19) and $η_{N} = \int η (ρ, Y_{e}, T, ε) ε^{- 1} d ε .$ ${\eta _N} = \int \eta \left( {\rho ,{Y_e},T,\varepsilon } \right){\varepsilon ^{ - 1}}{\rm{d}}\varepsilon .$ (22)

This definition ensures that neutrinos are in thermal equilibrium with the fluid when the optical depth (due to absorption) is high.

The energy-dependent emissivity, absorption opacity, and scattering opacity used to obtain their gray counterparts were calculated using NuLib (O’Connor 2015; Sullivan et al. 2016). A Python script to generate energy-integrated NuLib tables is available in the NuLib git repository².

3.4 Neutrino flux energy

To evaluate the number-density flux for the closure of Eq. (4), it is necessary to estimate the flux-weighted average energy in Eq. (8). The flux-weighted average energy is more or less equal to the average energy of the neutrinos in optically thin regions. However, the flux-weighted average energy and the average energy can differ significantly in optically thick regions. The neutrino opacities increase rapidly with energy, and the optical depth is, therefore, lower for low-energy neutrinos. The effect of this energy dependence in the optical depth is that low-energy neutrinos can escape while high-energy neutrinos remain trapped. In other words, the average energy of the neutrino flux can be significantly lower than the average neutrino energy. In regions with high absorption opacity, neutrinos remain in equilibrium with the fluid, and the exact value of flux-weighted average energy is relatively unimportant. On the other hand, the average energy of the neutrino flux in Eq. (8) plays an important role in regions where the scattering opacity is high and the absorption opacity is low (Foucart et al. 2016).

Due to the nature of a gray scheme, an energy-dependent scattering effect can not be fully captured with our approach, but simply ignoring it leads to a severe overestimate of the neutrino energy and underestimates the neutrino number diffusion rate. To obtain an estimate for the flux-weighted average energy Foucart et al. (2016) started from the fact that in the purely diffusive limit, assuming that the opacities go as ε², one expects (Rosswog & Liebendoerfer 2003) $〈 ε^{F} 〉 = \frac{F_{1} (η_{v})}{F_{0} (η_{v})} T_{v} .$ $\left\langle {{\varepsilon ^F}} \right\rangle = {{{F_1}\left( {{\eta _v}} \right)} \over {{F_0}\left( {{\eta _v}} \right)}}{T_v}.$ (23)

Based on the expected flux-weighted average energy, Foucart et al. (2016) introduced two new scalars and wrote down a simple equation for modeling the ratio of the flux-weighted average energy and the average neutrino energy: $\frac{〈 ε^{F} 〉}{〈 ε 〉} = \frac{F_{3} (η_{v}) F_{0} (η_{v}) - s^{F} (F_{3} (η_{v}) F_{0} (η_{v}) - F_{2} (η_{v}) F_{1} (η_{v}))}{F_{3} (η_{v}) F_{0} (η_{v}) - s^{C} (F_{3} (η_{v}) F_{0} (η_{v}) - F_{2} (η_{v}) F_{1} (η_{v}))} .$ ${{\left\langle {{\varepsilon ^F}} \right\rangle } \over {\langle \varepsilon \rangle }} = {{{F_3}\left( {{\eta _v}} \right){F_0}\left( {{\eta _v}} \right) - {s^F}\left( {{F_3}\left( {{\eta _v}} \right){F_0}\left( {{\eta _v}} \right) - {F_2}\left( {{\eta _v}} \right){F_1}\left( {{\eta _v}} \right)} \right)} \over {{F_3}\left( {{\eta _v}} \right){F_0}\left( {{\eta _v}} \right) - {s^C}\left( {{F_3}\left( {{\eta _v}} \right){F_0}\left( {{\eta _v}} \right) - {F_2}\left( {{\eta _v}} \right){F_1}\left( {{\eta _v}} \right)} \right)}}.$ (24)

Here s^C represents the fraction of neutrinos that have passed through a significant optical depth. s^F is a scalar that allows the average energy of the neutrino flux to be reduced in regions where s^C << 1. To ensure that s^C goes toward zero in the optically thick region, where the neutrinos are in equilibrium with the fluid, and toward s^F in the optically thin regions, where the neutrinos are streaming, one possible choice is $s^{C} = \frac{N s_{0}^{C} + α ℱ s^{F} d t}{N + α ℱ d t + {\bar{η}}_{N} α d t} .$ ${s^C} = {{Ns_0^C + \alpha {\cal F}{s^F}{\rm{d}}t} \over {N + \alpha {\cal F}{\rm{d}}t + {{\bar \eta }_N}\alpha {\rm{d}}t}}.$ (25)

In Eq. (25), s₀^C is the value at the beginning of the time step and 𝓕 is an estimate for the neutrino number flux. The scalar s_F can be defined as follows $s^{F} = \frac{s^{C} + τ}{1 + τ} .$ ${s^F} = {{{s^C} + \tau } \over {1 + \tau }}.$ (26)

The optical depth can be estimated based on the flux factor $(ξ = \sqrt{H^{μ} H_{μ}} / J)$ $(\xi = \left. {\sqrt {{H^\mu }{H_\mu }} /J} \right)$ and a free parameter β as follows: $ξ = \frac{1}{1 + β τ} .$ $\xi = {1 \over {1 + \beta \tau }}.$ (27)

The approximate neutrino number flux was set to $ℱ = ξ N {(\frac{F_{3} (η_{v}) F_{0} (η_{v}) - s^{F} (F_{3} (η_{v}) F_{0} (η_{v}) - F_{2} (η_{v}) F_{1} (η_{v}))}{F_{2} (η_{v}) F_{1} (η_{v})})}^{2} .$ ${\cal F} = \xi N{\left( {{{{F_3}\left( {{\eta _v}} \right){F_0}\left( {{\eta _v}} \right) - {s^F}\left( {{F_3}\left( {{\eta _v}} \right){F_0}\left( {{\eta _v}} \right) - {F_2}\left( {{\eta _v}} \right){F_1}\left( {{\eta _v}} \right)} \right)} \over {{F_2}\left( {{\eta _v}} \right){F_1}\left( {{\eta _v}} \right)}}} \right)^2}.$ (28)

4 Comparison of Foucart et al. (2016) to energy-dependent transport

To evaluate the accuracy of the original gray scheme, we performed a set of tests in one dimension (1D). The tests were all performed with β = 6 (in Eq. (27)). The choice of β = 6 was based on the recommendations of Foucart et al. (2016) and initial tests which indicated that this value lead to reasonable results. We started by considering the post-bounce evolution of a 30 solar-mass progenitor (Sukhbold et al. 2018) using the SFHo equation of state (EOS; Steiner et al. 2013). The progenitor is first evolved until 20 ms post bounce, using energy-dependent neutrino transport. From that point on, the fluid velocity is set to zero, and the rest of the fluid properties are assumed to be constant in time. Under these conditions, the system settles into a steady state. The results of the hydrostatic simulations are summarized in Table 1, where the first line represents the results from the simulation with fully energy-dependent neutrino transport, the second line shows the results when using the prescription of Foucart et al. (2016), and the last row shows the results of our improved gray transport (which we return to below). The original gray transport shows relatively good agreement with the energy-dependent neutrino transport. The average energies of all neutrino species agree to within ~20% in the two methods. The luminosities of the electron-type neutrinos agree within ~10%, the electron anti-neutrino luminosities agree to ~20%, and the heavy-lepton neutrino luminosities differ by almost 40%. The luminosities in the simulation using the original gray scheme are systematically smaller than those from the simulation with the energy-dependent transport.

O’Connor et al. (2018) compared results from several core-collapse supernova simulation codes. The initial conditions and input physics were carefully controlled to be as equal as possible across all the codes. Comparing 1D simulations of the first 500ms post bounce of a 20 solar mass progenitor, O’Connor et al. (2018) reported good agreement in key diagnostic quantities for all the codes. In Figs. 1–3 we compare our result for the same 20 solar mass progenitor with the results of O’Connor et al. (2018). The results presented in O’Connor et al. (2018) are represented by semitransparent colored lines. The solid purple line shows the results from the simulation with our implementation of the gray transport from Foucart et al. (2016). The luminosities in the simulation with the transport from Foucart et al. (2016) are initially significantly lower than the other simulations but agree better at later times (see Fig. 1). For the electron neutrinos and electron anti-neutrinos, the luminosities remain consistently lower in our test simulation than in the rest of the simulations. The luminosity of the heavy-lepton neutrinos eventually overtakes the luminosities reported in O’Connor et al. (2018) and remains high toward the end of the simulation. The transport of Foucart et al. (2016) leads to average neutrino energies for the electron neutrinos and electron anti-neutrinos that are within ~20% of the values in the rest of the simulations (see the left panel of Fig. 2). Additionally, we see that the average energy of the heavy-lepton neutrinos is ~60% larger compared to the typical value seen in the other simulations (see the right panel of Fig. 2). From Fig. 3, we see that both the PNS radius and the shock radius are significantly larger in the gray simulations compared to the values from the other simulations. The scheme proposed by Foucart et al. (2016) leads to a PNS which is typically 10 km to 20 km larger and a shock radius that is more than 30 km larger than what we expect based on the results of O’Connor et al. (2018).

The 1D tests illustrate that the scheme by Foucart et al. (2016) produces reasonable results, but differs significantly from the energy-dependent transport in terms of PNS properties, shock evolution, and neutrino emission. Consequently, we observe qualitative differences in the dynamical evolution of the 1D simulations with energy-dependent transport and our implementation of the original gray neutrino transport.

Table 1

Neutrino luminosities and average neutrino energy in hydrostatic simulations of the post bounce phase of a 30 solar mass progenitor.

5 Improved gray neutrino transport

5.1 Updated neutrino opacities

The choice of energy-averaged neutrino opacities is a natural place to start when investigating the root cause of the discrepancies we observed in our tests. There are two apparent potential issues with the opacities defined by Eq. (19). The accuracy of extrapolating the opacities based on the temperature ratio will decrease as the difference between the neutrino temperature and the fluid temperature increases. Furthermore, Eq. (18) does not account for the fact that the term κ_H in the evolution equations should be averaged with the neutrino momentum density (see Eq. (17)).

Since the neutrino transport and the hydrodynamics of any given simulations are inherently nonlinear, it can be challenging to construct comparisons that isolate a single component in dynamical simulations. We, therefore, based our initial evaluation of the accuracy of the correction applied to the opacities on the hydrostatic 1D simulations. The properties of the neutrino field, in terms of the neutrino temperature, enters into Eq. (19). Consequently, a direct comparison between simulations with energy-dependent and gray neutrino transport is complicated, since any change to the neutrino opacities leads to a change in the system’s steady state.

To isolate the test of the averaging procedure we performed the following evaluation: Using the spectral properties of the radiation fields in the simulation with energy-dependent neutrino transport, we computed $κ_{a}^{ed} = \frac{\int J (ε) κ_{a} (ρ, Y_{e}, T, ε) d ε}{\int J (ε) d ε},$ $\kappa _a^{{\rm{ed }}} = {{\int J (\varepsilon ){\kappa _a}\left( {\rho ,{Y_e},T,\varepsilon } \right){\rm{d}}\varepsilon } \over {\int J (\varepsilon ){\rm{d}}\varepsilon }},$ (29)

and $κ_{H}^{ed} = \frac{\int [κ_{a} (ρ, Y_{e}, T, ε) + κ_{s} (ρ, Y_{e}, T, ε)] H^{r} (ε) d ε}{\int H^{r} (ε) d ε},$ $\kappa _H^{{\rm{ed }}} = {{\int {\left[ {{\kappa _a}\left( {\rho ,{Y_e},T,\varepsilon } \right) + {\kappa _s}\left( {\rho ,{Y_e},T,\varepsilon } \right)} \right]} {H^r}(\varepsilon ){\rm{d}}\varepsilon } \over {\int {{H^r}} (\varepsilon ){\rm{d}}\varepsilon }},$ (30)

which we compared to the corresponding opacities calculated from Eq. (19). We obtained the neutrino temperature in Eq. (19) from the energy-dependent simulation by fitting the neutrino spectrum to a blackbody function. The fitting procedure inevitably induces errors in the neutrino temperature, but it is more consistent than comparing two different steady states.

In Figs. 4 and 5, we show the opacity values for κ_a and κ_H, respectively. In both figures, the red lines represent the calculations with the full neutrino spectrum (as per Eqs. (29) and (30)), and the green lines indicate the results from the equilibrium assumption (following Eq. (19)). The solid blue lines represent the new way of calculating the opacities we developed in this work, which we return to later in this section. The figures show κ_a and κ_H for the three different neutrino species as a function of radius. The top panel represents the electron neutrinos, the middle the electron anti-neutrinos, and the bottom panel shows the heavy-lepton neutrino species. The opacities for the heavy-lepton neutrinos refer to effective opacities where the pair-production processes have been included in a simplified manner (see Sect. 2).

At low radii, when the density is high and the neutrinos are trapped, we find good agreement between the average absorption opacities calculated directly from the neutrino spectrum (red curves) and gray opacities with the correction from Foucart et al. (2016; green curves). We see discrepancies as the density drops and the neutrinos start to decouple. The difference between the absorption for the heavy-lepton neutrinos is particularly large, see bottom panel of Fig. 4. In Fig. 5, the green line representing the blackbody assumption consistently yields higher values of κ_H than the red line based on the full neutrino spectrum.

The correction factor in Eq. (19) is an extrapolation away from equilibrium between the fluid and neutrinos. The anchor point $κ_{a / H}^{e q}$ $\kappa _{a/H}^{eq}$ is calculated assuming that the neutrino temperature is equal to the fluid temperature. This approach is reasonable as long as the neutrino temperature is close to the fluid temperature, but we expect a decrease in accuracy once the neutrinos decouple from the fluid. To improve the precision of the opacities, we used an interpolation method with pre-computed values of the opacities defined as follows $κ_{a / H} (T_{v}) = \frac{\int B_{v} (T_{ν}, μ_{v}, ε) κ_{a / H} (ρ, Y_{e}, T, ε) d ε}{\int B_{v} (T_{v}, μ_{v}, ε) d ε} .$ ${\kappa _{a/H}}\left( {{T_v}} \right) = {{\int {{B_v}} \left( {{T_\nu },{\mu _v},\varepsilon } \right){\kappa _{a/H}}\left( {\rho ,{Y_e},T,\varepsilon } \right){\rm{d}}\varepsilon } \over {\int {{B_v}} \left( {{T_v},{\mu _v},\varepsilon } \right){\rm{d}}\varepsilon }}.$ (31)

We chose a set of 12 linearly spaced neutrino temperatures from 0.1 MeV to l6 MeV. The neutrino opacities for each of the 12 neutrino temperatures were tabulated and stored in a four-dimensional table. We linearly interpolated from the tabulated values during the simulations. Thus in effect, the neutrino-interaction table remains the same size as in the energy-dependent simulations. Neutrino temperatures in the core can far exceed 16 MeV, but at that point the neutrinos are trapped and in equilibrium with the fluid. We found that reverting to Eq. (19) for the opacities when T_v > 16 MeV was beneficial from a numerical standpoint.

The interpolation method show a notable improvement over Eq. (19) for heavy-lepton neutrinos. On the other hand, for electron neutrinos and electron anti-neutrinos, the results are qualitatively similar to the approach prescribed by Foucart et al. (2016). We also find that the interpolation approach solves an issue we had noticed in our implementation of the original gray scheme. We found an artificial increase in the Y_e of the material in front of the shock. The increase in Y_e was due to the large difference in temperature between the neutrinos and the fluid ahead of the shock, which caused an overestimation of the absorption rates ahead of the shock. Interpolating, instead of extrapolating, the opacities removed the unphysical increase in Y_e in the in-falling material.

After implementing the interpolation for the neutrino opacities, a different problem became apparent. Reducing the absorption opacity of the heavy-lepton neutrinos moves the neutri-nosphere deeper into the PNS, increasing the average energy of the emitted neutrinos. After decoupling, in terms of absorption, the heavy-lepton neutrinos still have to transverse a region of significant scattering opacities. The higher average energy increases the scattering and prolongs the time it takes the heavy-lepton neutrinos to escape the optically semi-thick region in the outer layers of the PNS, which in turn leads to larger PNS radii.

Since the opacities decrease rapidly with energy, the average energy of the flux is lower in the optically thick regions (low-energy neutrinos escape while high-energy neutrinos remain trapped). Consequently, averaging with H^r instead of J leads to smaller energy-averaged opacities. The effect of this is demonstrated in Fig. 5. In the PNS, the energy-averaged opacities are roughly reduced by a factor of 0.6–0.8 when averaging with H^r rather than with J, compare the red and green lines³. Since we do not have information regarding the energy dependence of the momentum flux in our energy-integrated scheme, we multiply κ_H with 0.6 in the core to mimic this effect.

Implementing the correction to κ_H led to improved agreement between the PNS and neutrino properties of the gray scheme and the full energy-dependent neutrino transport. However, the gray models overestimated the heating rate. We account for this by simply multiplying the opacities by 0.95 in the post-shock layer. We rely on the flux factor to transition between the two regimes. Since the flux factor lies between 0 and 1, it is convenient to define κ_α and κ_H as follows: ${\bar{κ}}_{a} = (\cos^{2} (\frac{π ξ}{2}) + 0.95 {sin}^{2} (\frac{π ξ}{2})) κ_{a} (T_{v}),$ ${\bar \kappa _a} = \left( {{{\cos }^2}\left( {{{\pi \xi } \over 2}} \right) + 0.95{{\sin }^2}\left( {{{\pi \xi } \over 2}} \right)} \right){\kappa _a}\left( {{T_v}} \right),$ (32)

and ${\bar{κ}}_{H} = (0.6 {cos}^{2} (\frac{π ξ}{2}) + 0.95 {sin}^{2} (\frac{π ξ}{2})) κ_{H} (T_{v}) .$ ${\bar \kappa _H} = \left( {0.6{{\cos }^2}\left( {{{\pi \xi } \over 2}} \right) + 0.95{{\sin }^2}\left( {{{\pi \xi } \over 2}} \right)} \right){\kappa _H}\left( {{T_v}} \right).$ (33)

Here κ_a(T_V) and κ_H(T_v) refers to the interpolated values of the κ_a and κ_H, respectively. We note that the factor of 0.6 does not appear in absorption opacity, which shows up in the term W[η − κ_aJ] in the evolution equations, the correction factor should go to one in the optically thick limit to ensure detailed balance. The exact choice of the functional form of ${\bar{κ}}_{a} and {\bar{κ}}_{H}$ ${\bar \kappa _a}{\rm{ and }}{\bar \kappa _H}$ is somewhat arbitrary and several other good choices are likely to exist. We tested a few different variations before settling on the Eq. (32) and Eq. (33).

The opacities defined in Eqs. (32) and (33), represented by the blue lines in Figs. 4 and 5, respectively. For the heavy-lepton neutrinos, we see good agreement between ${\bar{κ}}_{a} and κ_{a}^{ed}$ ${\bar \kappa _a}{\rm{ and }}\kappa _a^{{\rm{ed}}}$ . The absorption opacities for the electron neutrinos and electron anti-neutrinos agree relatively well between all three methods, but ${\bar{κ}}_{a}$ ${\bar \kappa _a}$ is slightly closer to $κ_{a}^{ed} than κ_{a}^{eq}$ $\kappa _a^{{\rm{ed}}}{\rm{ than }}\kappa _a^{{\rm{eq}}}$ . For κ_H, we see that ${\bar{κ}}_{H}$ ${\bar \kappa _H}$ agrees better with the value from the energy-dependent calculation than $κ_{H}^{eq}$ $\kappa _H^{{\rm{eq}}}$ . However, there is still a some discrepancy between ${\bar{κ}}_{H} and κ_{H}^{ed}$ ${\bar \kappa _H}{\rm{ and }}\kappa _H^{{\rm{ed}}}$ , for all neutrino species.

Our tests showed that getting better results for κ_H is possible by choosing a different function and different coefficients in Eq. (33). However, the results shown in Figs. 4 and 5 are for one hydrostatic simulation. We found that tuning the scheme to fit one specific case led to worse results in other tests. In the end, Eqs. (32) and (33) were chosen to produce the best overall agreement in all our tests (see Sects. 5.3, and 6).

Fig. 1

Neutrino luminosities for the gray simulation compared to the simulations presented in O’Connor et al. (2018; semitransparent colored lines). The purple lines indicate the results from the simulation based on the original gray transport from Foucart et al. (2016), and the black lines represent the result from the simulation using the improvements detailed in this work. The luminosity of the electron neutrinos (solid lines) and electron anti-neutrinos (dot-dashed lines) are shown in the left panel, and the luminosity of the heavy-lepton neutrinos is shown in the right panel (dashed lines). Time is given in seconds after bounce. The inset in the left panel shows a zoom-in of the peak centered around ~0.l s after bounce.

Fig. 2

Average neutrino energies for the gray simulation (bright lines) compared to the simulations presented in O’Connor et al. (2018; semitransparent colored lines). The purple lines indicate the results from the simulation based on the original gray transport from Foucart et al. (2016), and the black lines represent the result from the simulation using the improvements detailed in this work. The average energy of the electron neutrinos (solid lines) and electron anti-neutrinos (dot-dashed lines) are shown in the left panel, and the average energies of the heavy-lepton neutrinos are shown in the right panel (dashed lines). Time is given in seconds after bounce.

Fig. 3

Shock radius (solid lines) and PNS radius (dashed lines) for the gray simulation (purple lines) compared to the simulations presented in O’Connor et al. (2018; semitransparent colored lines). The purple lines indicate the results from the simulation based on the original gray transport from Foucart et al. (2016) and the black lines represent the result from the simulation using the improvements detailed in this work. Time is given in seconds after bounce. The PNS radius is defined to be where the density drops below 10¹¹ g cm⁻³.

Fig. 4

Energy-averaged absorption opacity as a function of radius for a hydrostatic evolution of a 30 solar mass progenitor using the SFHo EOS. The dashed red lines show the opacities calculated with the neutrino spectrum from the energy-dependent neutrino transport (Eq. (29)). The dashed green lines show the opacities calculated according to Eq. (19). The solid blue lines represent the value according to Eq. (32). The panels show results for electron neutrinos (top), electron anti-neutrinos (middle), and heavy-lepton neutrinos (bottom). The opacities for the heavy-lepton neutrinos refer to effective opacities, which incorporate pair-production processes in a simplified way.

Fig. 5

Sum of the energy-averaged absorption and energy-averaged scattering opacities as a function of radius for a hydrostatic simulation of a 30 solar mass progenitor using the SFHo EOS. The dashed red lines show the opacities calculated with the neutrino spectrum from the energy-dependent neutrino transport (Eq. (30)). The dashed green lines show the opacities calculated according to Eq. (19). The blue lines represent the value according to Eq. (33). The panels show results for electron neutrinos (top), electron anti-neutrinos (middle), and heavy-lepton neutrinos (bottom). The opacities for the heavy-lepton neutrinos refer to effective opacities, which incorporate pair-production processes in a simplified way.

5.2 Updated flux-weighted average neutrino energy

To improve the machinery described in Sect. 3.4, we begin with the initial assumption of the ratio of the average neutrino energy and the flux-weighted average neutrino energy given in Eq. (23). The assumption is valid in the diffusion limit, if the opacities go as energy squared. In the free-streaming limit we would expect the energy of the flux to be the same as the average energy of the neutrino field. In reality, the radiation is neither purely diffusive nor purely free-emission, which means that $\frac{F_{1} (η_{v})}{F_{0} (η_{v})} T_{v} < 〈 ε_{F} 〉 < \frac{F_{3} (η_{v})}{F_{2} (η_{v})} T_{v} .$ ${{{F_1}\left( {{\eta _v}} \right)} \over {{F_0}\left( {{\eta _v}} \right)}}{T_v} < \left\langle {{\varepsilon _F}} \right\rangle < {{{F_3}\left( {{\eta _v}} \right)} \over {{F_2}\left( {{\eta _v}} \right)}}{T_v}.$ (34)

Results from the energy-dependent transport indicate that the flux-weighted average energy matches well with $〈 ε_{F} 〉 = \frac{F_{2} (η_{v})}{F_{1} (η_{v})} T_{v},$ $\left\langle {{\varepsilon _F}} \right\rangle = {{{F_2}\left( {{\eta _v}} \right)} \over {{F_1}\left( {{\eta _v}} \right)}}{T_v},$ (35)

especially for the heavy-lepton neutrinos (in the region with high scattering opacity) for which the second term in Eq. (8) is particularly important. We, therefore, change the order of the Fermi integrals in Eqs. (24) and (28) to better match the results from the energy-dependent transport. The updated expressions are $ℱ = ξ N {(\frac{F_{3} (η_{v}) F_{1} (η_{v}) - s^{F} (F_{3} (η_{v}) F_{1} (η_{v}) - F_{2} (η_{v}) F_{2} (η_{v}))}{F_{2} (η_{v}) F_{2} (η_{v})})}^{2},$ ${\cal F} = \xi N{\left( {{{{F_3}\left( {{\eta _v}} \right){F_1}\left( {{\eta _v}} \right) - {s^F}\left( {{F_3}\left( {{\eta _v}} \right){F_1}\left( {{\eta _v}} \right) - {F_2}\left( {{\eta _v}} \right){F_2}\left( {{\eta _v}} \right)} \right)} \over {{F_2}\left( {{\eta _v}} \right){F_2}\left( {{\eta _v}} \right)}}} \right)^2},$ (36)

and $\frac{〈 ε^{F} 〉}{〈 ε 〉} = \frac{F_{3} (η_{v}) F_{1} (η_{v}) - s^{F} (F_{3} (η_{v}) F_{1} (η_{v}) - F_{2} (η_{v}) F_{2} (η_{v}))}{F_{3} (η_{v}) F_{1} (η_{v}) - s^{C} (F_{3} (η_{v}) F_{1} (η_{v}) - F_{2} (η_{v}) F_{2} (η_{v}))} .$ ${{\left\langle {{\varepsilon ^F}} \right\rangle } \over {\langle \varepsilon \rangle }} = {{{F_3}\left( {{\eta _v}} \right){F_1}\left( {{\eta _v}} \right) - {s^F}\left( {{F_3}\left( {{\eta _v}} \right){F_1}\left( {{\eta _v}} \right) - {F_2}\left( {{\eta _v}} \right){F_2}\left( {{\eta _v}} \right)} \right)} \over {{F_3}\left( {{\eta _v}} \right){F_1}\left( {{\eta _v}} \right) - {s^C}\left( {{F_3}\left( {{\eta _v}} \right){F_1}\left( {{\eta _v}} \right) - {F_2}\left( {{\eta _v}} \right){F_2}\left( {{\eta _v}} \right)} \right)}}.$ (37)

The scalar s^C represents the fraction of neutrinos that have gone through a significant optical depth since being emitted. This ratio reaches values close to one well before the emitted neutrinos have escaped the PNS. Since s^C = 1 implies that 〈ε^F〉 = 〈ε〉 (see Eq. (24)), the energy correction is effectively turned off once s^C reaches values close to unity. The fact that s^C becomes one before the neutrinos have escaped the scattering region leads to an underestimate of the reduction of the average energy of the flux. We, therefore, set $s^{C} = \frac{N s_{0}^{C} + α ℱ s^{F} d t}{N (1 + τ) + α ℱ d t + {\bar{η}}_{N} α d t} .$ ${s^C} = {{Ns_0^C + \alpha {\cal F}{s^F}{\rm{d}}t} \over {N(1 + \tau ) + \alpha {\cal F}{\rm{d}}t + {{\bar \eta }_N}\alpha {\rm{d}}t}}.$ (38)

The only difference between Eqs. (38) and (25) is the term N(1 + τ) in the denominator, which serves to take into account how much optical depth the neutrinos have yet to travel through. Effectively, our choice of s^C goes to one much slower than Eq. (25). Furthermore, we found that the estimate of the optical depth suggested by Foucart et al. (2016) underestimates the optical depth. We, therefore, calculate τ as follows $τ (r) = \int_{r}^{\infty} {〈 {\bar{κ}}_{H} 〉}_{Ω}^{- 1} d r,$ $\tau (r) = \int_r^\infty {\left\langle {{{\bar \kappa }_H}} \right\rangle _\Omega ^{ - 1}} {\rm{d}}r,$ (39)

where ${〈 {\bar{κ}}_{H} 〉}_{Ω}$ ${\left\langle {{{\bar \kappa }_H}} \right\rangle _\Omega }$ represents the angular average of ${\bar{κ}}_{H}$ ${\bar \kappa _H}$ at any given radius. We use ${\bar{κ}}_{H}$ ${\bar \kappa _H}$ to calculate the optical depth since we are attempting to estimate the average energy of the neutrino flux. Directly calculating τ is advantageous since it completely does away with the β parameter of Foucart et al. (2016), a parameter that influences the original scheme’s results. Foucart et al. (2016) reported a 5–10% variation in the neutrino luminosites and average energies for reasonable values of β (∈ (0, 8]). We found a similar dependence on the β parameter, as Foucart et al. (2016), in our implementation of their scheme. While calculating the optical depth directly improves the result of our gray scheme, relying on the approximate value from Eq. (27) is not detrimental to the accuracy of the transport. Furthermore, Eq. (39) assumes that τ is well represented by a spherically symmetric radial profile. Core-collapse supernovae are fairly close to spherical, but the implicit assumption the Eq. (39) is likely not appropriate for binary mergers.

The accuracy of our gray scheme depends on the two free parameters we introduced in Eqs. (32) and (33). Consequently, it is important to chose appropriate values for the two parameters. In our tests, we found that the parameter that appears in the first term of Eq. (33) should lie between 0.5 and 0.8, changes of 0.2 results in noticeably different PNS properties. The second parameter is important for the outcome of the supernova simulations, explosion versus non-explosion, and is sensitive to changes on the order of 5%.

5.3 Numerical tests

After implementing the changes detailed above, we performed the same 1D tests we initially used to test the scheme of Foucart et al. (2016). The last line in Table 1 represents the results of our updated scheme. Except for the average energy of the heavy-lepton neutrinos, we find agreement with the energy-dependent scheme within 5–10%. The average energy of the heavy-lepton neutrinos agrees within 30%. Since we cannot capture the effects of inelastic scattering and struggle to model the full impact of scattering screens, we do not expect to be able to reproduce the average energy of the heavy-lepton neutrinos well with a gray scheme.

The black curves in Figs. 1–3 show the results using the updated gray scheme. The new gray scheme generally leads to better agreement with the other simulations in all quantities except the average energy of the heavy-lepton neutrinos. The average shock radius is still larger than the typical value seen in the original simulations, but the PNS radius agrees well with the results reported in O’Connor et al. (2018). Getting the properties of the PNS right is important for predicting the GWs emitted by core-collapse supernovae (Müller et al. 2013; Sotani & Takiwaki 2016; Sotani et al. 2017, 2021; Andresen et al. 2019; Morozova et al. 2018; Westernacher-Schneider et al. 2019; Eggenberger Andersen et al. 2021; Vartanyan et al. 2022; Mori et al. 2023).

6 Two-dimensional tests

We extended our testing of the gray scheme by conducting 2D axisymmetric core-collapse supernova simulations, using both the fully energy-dependent neutrino transport and the gray transport. The energy-dependent transport was carried out using 12 energy groups, inelastic scattering, and energy-bin coupling effects. The simulations are grouped into two sets. The first set is composed of five progenitors with ZAMS masses of 11, 12, 13, 14, and 15 solar masses (Sukhbold et al. 2018). We simulated the post-bounce evolution of each progenitor three times, one for each neutrino transport method. We used the SFHo EOS, and the finest grid resolution was 488 m. The progenitors were evolved until 20 ms after bounce in 1D with energy-dependent neutrino transport. The 2D simulations were initiated from the post bounce 1D profiles. The simulations are labeled with an “s” followed by the progenitor ZAMS mass. For each progenitor, the average shock radius and the PNS radius in the corresponding simulations are shown in Fig. 6.

The second set of simulations is based on a single progenitor with a ZAMS mass of 20 solar masses (Woosley & Heger 2007). We conducted the simulations using five different versions of the SRO EOS (Schneider et al. 2017, 2019). The five EOS used are chosen from a larger set (of 96) that span a parameter space allowed by terrestrial experiments and astrophysical observations. Specifically, the five different versions of the EOS lead to different simulation outcomes in energy-dependent simulations, which enables us to check if our gray scheme reproduces this result and to test different regimes. The finest grid resolution was 325 m. The energy-dependent simulations from the second set have already been published in Eggenberger Andersen et al. (2021); we refer to their work for a detailed analysis of the simulations. For each EOS, the average shock radius and the PNS radius in the corresponding simulations are shown in Fig. 7. We observe a notable offset in the initial PNS radius of the simulations with gray transport compared to the energy-dependent transport simulations. The offset is related to the fact that the original simulations were evolved through bounce in 2D, while the gray simulations were mapped from a spherically symmetric profile 20 ms after bounce. It is necessary to start the gray simulations post bounce, because energy dependent effects play an important role during the collapse phase. The absence of prompt convection, caused by evolving the model through bounce in 1D, influences the structure of the PNS and leads to a smaller initial PNS radius than when evolving the collapse phase in 2D. The difference is transient, and the excellent agreement of the shock radii at early times demonstrates that this effect is relatively unimportant for the further dynamic evolution of the system.

In general, we find remarkably good qualitative agreement between simulations using our gray transport and simulations using energy-dependent neutrino transport for a wide range of progenitor and for several different EOSs (see Figs. 6 and 7). We present a side-by-side comparison of the simulations from the first set in Fig. 6. We plot the average shock radius and the PNS radius (we define the PNS surface as the radius where the density drops below 10¹¹ g cm⁻³) for each progenitor. The left column displays the shock radius, while the right column shows the PNS radius. Each row of Fig. 6 corresponds to a different progenitor. The blue solid line represents the gray scheme implemented in this work, the green dashed-dotted line represents the method by Foucart et al. (2016), and the red dashed line represents the simulations utilizing energy-dependent transport.

From the right-hand side panels of Fig. 6, we see that the PNS radii are systematically larger in the simulation using the prescription of Foucart et al. (2016). Consistent with the results from our 1D tests, the modifications we suggest lead to improved agreement between the energy-dependent and energy-integrated transport in terms of the PNS radii. However, the gray transport still leads to larger PNS radii at late times than the energy-dependent transport. For the s11 progenitor, at 400 ms after bounce the PNS radius is ~33 km, ~35.5 km, and ~38 km in the simulations using energy-dependent transport, our gray implementation, and the gray transport according to Foucart et al. (2016), respectively. For the s14 progenitor, again at 400 ms post bounce, we find a PNS radius of ~34 km using energy-dependent transport, ~36.5 km using our gray transport and ~39 km using the original gray transport. The rest of the progenitors show similar trends. As in 1D, the reduced PNS radii we observe after implementing our changes result from the reduced scattering opacities for heavy-lepton neutrinos. Consequently, the heavy-lepton neutrinos are emitted deeper in the PNS. Therefore, the updated gray transport leads to larger average energies for the heavy-lepton neutrinos than the original implementation.

In terms of explodability and time of shock-revival, we find good qualitative agreement between the energy-dependent transport and our gray transport. The gray scheme is more conducive to shock revival, and shock expansion tends to set in sooner in simulations using the gray transport than in simulations using energy-dependent transport. The situation is notably reversed for the original gray scheme; the original gray transport tends to lead to later shock revival than the energy-dependent transport.

We show the heating rate as a function of time for models s11, s12, s13, s14, and s15 in Fig. 8. The total heating rate is defined as the volume integral of the change in the specific internal energy due to neutrino-fluid interactions, in regions where the change is positive. The spatial integral is limited to regions where the density is less than 3 × 10¹⁰ g cm⁻³ and the entropy is greater than 6 k_B per baryon (k_B denotes Boltzmann constant). In Fig. 8, blue curves indicate results from simulations with our gray scheme, green curves show the results from simulations using the neutrino transport of Foucart et al. (2016), and red curves correspond to results from simulations with energy-dependent neutrino transport. Two trends are immediately visible in Fig. 8: 1) the neutrino transport of Foucart et al. (2016) leads to lower heating rates than the fully energy-dependent transport, 2) our gray scheme injects too much energy into the fluid at early times, the heating rates are to high, with the notable exception of model s12. The shock radii of all three simulations of the s12 progenitor follow a very similar trajectory prior to 200 ms post bounce, see the second panel from the top in the left column of Fig. 6. The large heating rates at early times, and the larger heating rate in general, is a likely reason for the differences in explodability we observe between our gray transport and the energy-dependent transport. Especially since the evolution of the gray s12 simulation, where this excessive early-time heating is absent, agrees so well with its energy-dependent counterpart.

In our work, we have focused on achieving the best possible agreement in the properties of the neutrino radiation, the properties of the PNS, and the evolution of the shock. The aforementioned quantities are better tracers for the global behavior of the simulations than the total heating rate. The heating rate does not account for cooling or where the heating occurs. Nevertheless, the high heating rates of the gray scheme indicate that our scheme could be improved in the future.

We compare the simulations from the second set in Fig. 7. We plot the average shock radius and the PNS radius for each realization of the post-bounce evolution of the s20 progenitor. The left column displays the shock radius, while the right column shows the PNS radius. Each row of Fig. 7 corresponds to a distinct EOS. We follow the labeling scheme on Eggenberger Andersen et al. (2021; see Table 2 for an overview of the naming convection we used for the simulations), where the progenitor name is followed by the value chosen for the effective mass of neutrons at nuclear saturation density, for symmetric nuclear matter, in the EOS. A second number following the effective mass represents the value chosen for the isoscalar incompressibility modulus (in units of MeV per baryon), if no number is given then the isoscalar is set to the baseline value of 230 MeV. The solid blue line represents the gray scheme implemented in this work, the green dashed-dotted line represents the method by Foucart et al. (201 ), and the red dashed line represents the simulations utilizing energy-dependent transport.

Similar to the simulations in the first set, for the second set of simulations based on the s20 progenitor, we find that the original gray scheme consistently results in larger PNS radii than the energy-dependent transport. Furthermore, we see that the updated gray transport again leads to larger PNS radii than the energy-dependent transport, from approximately ~200 ms post bounce.

For the s20_0.55 models we observe good agreement in the shock evolution for all three transport methods, and none of the simulations explode within the first 800 ms post bounce (see top right panel of Fig. 7). The s20_0.95 models produce successful explosions and we observe the same trend we did in the first set of simulations; the original scheme by Foucart et al. (201 ) leads to later shock revival and our gray scheme leads to earlier explosions compared to the energy-dependent transport. However, shock revival sets in within 50 ms in the simulation using our gray transport and the simulation with energy-dependent transport. After an initial period of expansion, the shock begins to retreat around 100 ms post bounce in the simulation with energy-dependent transport (bottom right panel of Fig. 7). This period of shock retraction is observed to a lesser degree in the model with the original gray transport. It is absent in the simulation that incorporates the modified gray scheme. Similar behavior can be seen in the s20_0.75_200 simulations.

The s20_0.75 simulations are interesting, Eggenberger Andersen et al. (2021) found very similar shock evolution for all values of the incompressibility modulus, but our gray transport produces qualitatively different behavior for all three models. On the other hand, the transport of Foucart et al. (201 ) yields better agreement across the three EOS variations. The 2D tests show that the gray transport shifts the outcome toward explosion for progenitors in the boundary region between explodability and non-explodability. Consequently, our findings warrant caution in using gray transport to predict the outcome of individual progenitors precisely. However, our gray transport successfully reproduces the main dependence of explodability on the EOS parameters reported in Schneider et al. (2019); Yasin et al. (2020); Eggenberger Andersen et al. (2021), namely that increasing the effective mass facilitates shock revival.

Fig. 6

Average shock radius (left column) and the PNS radius (right column) for simulations based on five progenitors with ZAMS masses of 11, 12, 13, 14, and 15 solar masses. We performed three simulations for each progenitor, using different neutrino transport methods. The dashed red line corresponds to simulations with fully energy-dependent transport, the green lines indicates simulations using the gray scheme of Foucart et al. (2016), and the blue lines show results from simulations using the gray transport proposed in this work. Each row corresponds to a given progenitor; the progenitor model is indicated in the top left corner of each row. The simulations were performed with the SFHo EOS.

Table 2

Overview of the model names for our second set of 2D simulations.

Fig. 7

Average shock radius (left column) and the PNS radius (right column) for simulations based on one progenitor with a ZAMS mass of 20 solar masses. For each EOS we performed three simulations, one for each neutrino transport method. The dashed red line corresponds to simulations with fully energy-dependent transport, the green lines indicate simulations using the gray scheme of Foucart et al. (2016), and the blue lines show results from simulations using the gray transport proposed in this work. Each row corresponds to a given version of the SRO EOS; EOS parameters (see text) are indicated by the model name in the top left corner of each row.

Fig. 8

Heating rate as a function of time for models s11, s12, s13, s14, and s15. The blue curves indicate the results from simulations with our gray scheme, the green curves show the results from simulations using the neutrino transport of Foucart et al. (201 ), and the red curves correspond to the results from simulations with energy-dependent neutrino transport. Each row correspond to one model; the model names are shown in the bottom left corner of each panel.

7 Gravitational waves and neutrino emission

In this section, we present the neutrino luminosity, average energy, and GW signal for the energy-dependent and gray simulations based on the s12 progenitor. We do not show results for the gray implementation of Foucart et al. (201 ). Our decision to concentrate on s12 is based on two reasons. First, a detailed examination of all the simulations would result in an overly long and extensive analysis. Second, and more crucially, we have verified that the s12 simulations are representative of the totality of simulations performed in this work.

The GW signals are extracted from the hydrodynamic simulations by post-processing the output data using the quadrupole formula. In the transverse-traceless gauge, the GW tensor can be expressed in terms of two independent components, h₊ and h_×, but only one component is nonzero in axisymmetric simulations. Far away from the source, at a distance D, in the slow-motion limit, and assuming an observer located in the equatorial plane of the simulations, we have $h_{+} = \frac{3}{2} \frac{G}{c^{4} D} {\ddot{Q}}_{33} .$ ${h_ + } = {3 \over 2}{G \over {{c^4}D}}{\ddot Q_{33}}{\rm{. }}$ (40)

Here is the c the speed of light and G is Newton’s constant. ${\ddot{Q}}_{33}$ ${\ddot Q_{33}}$ denotes the second-order time derivatives of the zz-component of the quadrupole moment and is given by ${\ddot{Q}}_{33} = \frac{4}{3} \frac{d}{d t} [\int ρ ɀ v_{ɀ} d^{3} x] .$ ${\ddot Q_{33}} = {4 \over 3}{{\rm{d}} \over {{\rm{d}}t}}\left[ {\int \rho z{v_z}{{\rm{d}}^3}x} \right]{\rm{. }}$ (41)

In this form, one of the time derivatives of the original definition has been eliminated to avoid numerical problems associated with second-order derivatives (Oohara et al. 1997; Finn 1989; Blanchet et al. 1990). In Eq. (41), v_z is the Cartesian velocity, z is the Cartesian coordinate, and ρ is the local fluid density. We calculate the second time derivative numerically using numpy.gradient (Harris et al. 2020).

We compute spectrograms by applying short-time Fourier transforms (STFT) to h₊. The STFT is computed with a scipy.signal.stft using a Blackman window (Virtanen et al. 2020). We normalize the STFT and take the logarithm before plotting. The normalization is chosen so that the logarithmic value lies in the (−∞,0] range and is the same for every spectrogram we show. Before applying the STFT, we filter the signals using high-pass and low-pass filters, removing any part of the signal below 25 Hz and above 5000 Hz.

In Fig. 9, we present the GW amplitude and spectrograms for the two simulations. The top row displays the GW amplitude, while the bottom row shows the spectrograms. The left column of each row corresponds to the results from the energy-dependent neutrino transport simulation, and the right column showcases the findings from the gray neutrino transport simulation.

The time-dependent signals of both models (top row of Fig. 9) exhibit the stochastic variations typically associated with core-collapse supernova GWs. However, the typical amplitudes are ~15 cm in both simulations. The strong burst of GW emission visible in the simulation with energy-dependent transport around 50 ms post bounce is associated with prompt convection and is absent in the simulation with gray transport because it was initiated from a 1D profile post bounce. At first glance, the spectrograms from both models are very similar. However, there is a small shift toward lower frequencies in the simulation employing gray neutrino transport, compared to the simulation with energy-dependent transport. This frequency shift is consistent with the differences in PNS properties observed in simulations that utilize different neutrino transport methods. The gray transport leads to a less compact PNS at late times than the energy-dependent transport. It is, therefore, expected that the central frequency of the GW signal will be lower (Murphy et al. 2009; Mueller et al. 2013; Sotani & Takiwaki 2016; Andresen et al. 2017; Andresen 2017; Sotani et al. 2017, 2019, 2021; Torres-Forné et al. 2018, 2019b,a; Morozova et al. 2018; Radice et al. 2019; Eggenberger Andersen et al. 2021; Mezzacappa et al. 2023; Wolfe et al. 2023; Rodriguez et al. 2023). Additionally, the power gap is observed in both models. Interestingly, the location of the power gap (Morozova et al. 2018; Eggenberger Andersen et al. 2021) is shifted toward lower frequencies in the gray transport simulation.

Figure 10 shows the neutrino luminosities and average neutrino energies from the two simulations of the s12 progenitor. The top panel of the figure shows the luminosities, while the bottom panel shows the average energies. The neutrino luminosity and average neutrino energies are extracted at a radius of 500 km, we apply a Savitzky-Golay filter (Virtanen et al. 2020) to smooth the curves before plotting. In both panels, the red lines represent the results obtained from the simulation that employed energy-dependent neutrino transport. The blue lines correspond to the results from the gray neutrino transport simulation. Solid lines indicate electron neutrinos, dashed lines represent electron anti-neutrinos, and dotted-dashed lines correspond to heavy-lepton neutrinos.

For the electron neutrinos and electron anti-neutrinos, there is excellent agreement between the simulations in terms of both average energies and luminosities (the luminosity curves trace each other very well with only a small off-set, see the first panel of Fig. 10). However, the gray neutrino transport induces notable differences for the heavy-lepton neutrinos compared to the energy-dependent transport. Initially, the gray transport produces approximately 10 MeV higher energies, with around 30% lower luminosities than those observed in the energy-dependent transport simulation. At late times, the luminosity agreement improves to about 15–20%.

The ~ 10 MeV offset in average-energy of the neutrinos heavy-lepton, between the simulations with gray and energy-integrated transport, is present for all our models. Furthermore, all of our gray 2D models show the same 15–30% reduction in the heavy-lepton neutrino luminosity, compared to the simulations with energy-dependent transport, as the s12 models discussed in this section. For the s12 models, the energy-dependent and gray transport produces similar average energies and luminosities for the electron neutrinos and anti-neutrinos. The good agreement observed for the s12 models carries over to all our 2D models.

In Sect. 5.3, we speculated that the differences in the average energy of the heavy-lepton neutrinos were, at least in part, caused by the inclusion of inelastic scattering in the energy-dependent transport. To test the effects of inelastic scattering, we redid the energy-dependent s12 simulation without inelastic scattering. We see marginal differences. The largest is the aforementioned impact on the average energy where we saw an increase in the average energy of the heavy-lepton neutrinos of approximately 5 MeV. We also observed a slightly smaller PNS, and a smaller average shock radius. The luminosities of the heavy-lepton neutrinos increased by 5–10% when inelastic scattering was not included, the luminosity of the two other species stayed more or less the same. The model without inelastic scattering follows the same general evolution as the one with, but is over all less energetic than the simulation with inelastic scattering. The gray simulation matches the simulation without inelastic scattering better than the one with inelastic scattering.

The differences we observe in the neutrino emission are consistent with the results in 1D test simulations, see Sect. 5. Furthermore, we found that the differences in the neutrino emission in the 2D simulations using the original scheme and the energy-dependent transport were the same as in 1D (see Sect. 5). The original scheme proposed by Foucart et al. (2016) results in larger PNS radii compared to energy-dependent transport schemes. Consequently, by the end of the simulation, this leads to a reduction in the central frequency of the gravitational wave emission by approximately 300 Hz.

Fig. 9

GWs from two simulations of the s12 progenitor. The right column shows the signals from the simulation with energy-dependent neutrino transport and the left column shows the corresponding simulation performed with our gray neutrino transport. The top row shows the strain as a function of time and the bottom row shows the square of the short time Fourier transform calculated with scipy.signal.stft (Virtanen et al. 2020). The color scale is logarithmic, and the plots were normalized by a common factor.

Fig. 10

Neutrino luminosities (top panel) and average energies (bottom panel) from two simulations of the s12 progenitor. The red lines indicate the results from the simulation with energy-dependent neutrino transport. The blue lines show the results from the simulation performed with our gray neutrino transport. The solid lines indicate electron neutrinos, dashed lines show electron anti-neutrinos, and dot-dashed lines correspond to heavy-lepton neutrinos.

8 Performance

To evaluate the performance of the scheme, we conducted two separate 2D simulations. Each simulation was run for a total of 1000 time steps, with the initial resolution of the finest refinement level set at 325 m. In our tests, we found that the gray neutrino transport speeds up the overall code execution by a factor of four compared to the energy-dependent transport. When looking specifically at the neutrino transport, the speed increase is even more significant, surpassing a factor of five. The radiation transport takes approximately three times longer than the hydro-dynamic calculation. For comparison, in the energy-dependent scheme, a radiation step is ~15 times more expensive than a hydrodynamic step.

9 Summary and conclusion

In this work, we performed a comprehensive comparison of the energy-integrated neutrino transport scheme of Foucart et al. (2016) with fully energy-dependent neutrino transport. We found that the method outlined by Foucart et al. (2016) yields reasonable shock evolutions in our tests, but PNS radii that are 15–20% larger compared to those obtained in corresponding simulations using energy-dependent neutrino transport. We also found, in 2D test simulations, that our implementation of the original gray transport reduced the explodability of the models compared to the energy-dependent simulations.

Based on results from numerical tests, we implemented several modifications to the original gray scheme. We implemented the following changes to the approach proposed by Foucart et al. (2016):

1.
We updated the definition of energy-averaged neutrino opacities. Foucart et al. (2016) computed the scattering and absorption opacities by first calculating the opacities under the assumption that the neutrinos are in equilibrium with the fluid and then extrapolating away from equilibrium using the ratio (T_ν/T)². We assumed that the neutrinos obey a black-body spectrum (which was also assumed by Foucart et al. 2016), and computed the opacities for a set of predefined neutrino temperatures to create an interpolation grid in T_ν. After interpolating the opacities, we applied a set of correction factors. These correction factors, in short, amount to multiplying the opacities with 0.95 in the gain layer and reducing the opacities that appear in the evolution equations as (κ_a + κ_s)Hⁱ with 0.6 in optically thick regions.
2.
We adjusted the procedure for estimating the average energy of the neutrino flux. Specifically, we changed the order of the Fermi integrals in the equations, modeling the ratio of the average energy of the neutrino flux to the average neutrino energy. For consistency, we also changed the updated order of the Fermi integrals that appear in the estimate for the neutrino number flux. We modified the scalar s^C by introducing an Nτ term in the denominator. We directly computed the optical depth, eliminating reliance on the β parameter.

Intuitively, it is reasonable that computing the opacities from a grid of values is more accurate than extrapolating from a single point. While interpolating led to improvements across all neutrino species, the largest improvement was in the accuracy of the heavy-lepton neutrino opacities. Furthermore, in our initial implementation of the original gray transport we found an artificial increase in the Y_e of the material ahead of the shock. The source of this build-up was the extrapolation of the absorption opacity in a regime where the fluid temperature was very low (<1 MeV) and the neutrino temperature was high. Interpolating the opacities, instead of extrapolating, fixed this problem. Our approach reproduces the results of the energy-dependent transport fairly well; in addition, simulations utilizing the gray scheme are around four times faster than those using the energy-dependent neutrino transport.

Interpolating the opacities was not enough to achieve good agreement with the opacities calculated by averaging with the energy-density spectrum from simulations with energy-dependent transport. In the evolution equations there are terms which are proportional to (κ_s + κ_a)Hⁱ, which means that the opacities in these terms should be averaged with the neutrino momentum density instead of the neutrino energy density. Effectively, this reduces the opacities in the optically thick regions. We found that multiplying the opacities with 0.6 in the thick regions mimicked this effect reasonably well. We found that these corrections were essential for the heavy-lepton neutrinos. Additionally, we observed that the energy-averaged opacities lead to overheating in the gain layer, which we solved by reducing the opacities by 5% in this region.

The methodology of Foucart et al. (2016) for tracking the average energy of the neutrino flux presupposes that the average energy of neutrinos traveling through an optically thick screen will be reduced from $\frac{F_{3} (η_{v})}{F_{2} (η_{v})} T_{v} to \frac{F_{1} (η_{v})}{F_{0} (η_{v})} T_{v}$ ${{{F_3}\left( {{\eta _v}} \right)} \over {{F_2}\left( {{\eta _v}} \right)}}{T_v}{\rm{ to }}{{{F_1}\left( {{\eta _v}} \right)} \over {{F_0}\left( {{\eta _v}} \right)}}{T_v}$ . Empirically, we found that the energy reduction is closer to $\frac{F_{2} (η_{v})}{F_{1} (η_{v})} T_{v}$ ${{{F_2}\left( {{\eta _v}} \right)} \over {{F_1}\left( {{\eta _v}} \right)}}{T_v}$ . Consequently, we changed the order of the Fermi integrals involved in tracking the momentum-weighted average energy. As for the previous point, we found that this has the most significant effect on the propagation of heavy-lepton neutrinos.

Tracking the flux-weighted average energy entails calculating the optical depth of the neutrinos. Foucart et al. (2016) suggests an estimate that we found to underestimate the actual value calculated directly in our simulations. We therefore chose to calculate the optical depth directly. Relative to the other improvements, the effect of this correction is small and can likely be left out if calculating the optical depth is numerically challenging.

The changes detailed above were developed in an iterative procedure, where we implemented changes and tested their impact by performing 1D and 2D tests. After arriving at our final set of changes to the original transport, we performed several 2D test simulations of our gray scheme. Overall, we found that for a range of progenitors and for different EOSs the gray neutrino transport successfully reproduces key quantities from the more intricate energy-dependent transport. We find good qualitative and reasonable quantitative agreement between simulations with our gray transport and simulations with multi-group M1 transport. However, the gray scheme is slightly more conducive to shock revival, and we conclude that it will likely push simulations residing in the boundary region between successful and failed supernovae over the edge toward explosion.

We performed two sets of 2D test simulations. In the first set we tested the accuracy of the gray transport over a range of progenitors, spanning ZAMS masses from 11 to 15 solar masses, and used the SFHo EOS for every simulation. We found good qualitative and quantitative agreement for the evolution of the PNS and the shock radius for all models, except for one progenitor in which shock expansion was set approximately 100 ms earlier in the gray simulation. The scheme by Foucart et al. (2016) resulted in significantly delayed explosions, compared to the energy-dependent transport simulations, for two out of the five progenitors. The original scheme also consistently led to PNS radii that were ~10 km too large.

In the second set of 2D simulations, we kept the progenitor constant and studied the effects of varying the EOS. Recomputing the simulations performed by Eggenberger Andersen et al. (2021) with our gray scheme showed that it is slightly more conducive to explosions than the energy-dependent simulations. The particular model used for the comparison with Eggenberger Andersen et al. (2021) has been used by several authors that have reported successful and unsuccessful explosions (Melson et al. 2015a). A likely conclusion is that the gray scheme favors explosions in cases close to the boundary between failed and successful supernova explosions. Again, we found that the original gray transport hampers shock revival and leads to extended PNSs.

One motivation for implementing a gray scheme into FLASH is the ability to produce a larger and more diverse set of predictions for the multi-messenger signals from core-collapse supernovae than currently available in the literature. We verified that our gray transport produces GW and neutrino signals that agree well with those of the energy-dependent transport. In Sect. 7 we compared the predicted GW and neutrino signals for one particular model using the gray transport with the corresponding model with energy-dependent transport. The GW signals agreed well, but we observed a small shift toward lower frequencies in the gray model, which we attributed to the larger PNS radius observed in the gray model (a tendency that is systematic to the gray models). Furthermore, we found relatively good agreement in the neutrino emission of the two models. The largest discrepancy was in the average energy and luminosity of the heavy-lepton neutrinos, which, considering the limitations of a gray scheme, was expected.

Since we implemented the gray transport into a version of FLASH, which already has an energy-dependent M1 scheme, we were able to assess the consequences of every aspect of the gray transport. Such a detailed comparison is critical to understanding the limitations of commonly used gray schemes. The standard choice for the neutrino opacities can be improved to achieve better agreement with the energy-dependent neutrino transport. The corrections we suggest were empirically determined, but constructed to respect the expected behavior in the optically thick and thin limits. Other choices for the neutrino opacities will likely work as well as our suggested corrections. Regardless of which choice one makes, our recommendation is to validate any given prescription for the energy-averaged opacities, because strong nonlinear feedback couples the entire problem.

Acknowledgements

We thank Francois Foucart and Federico Schianchi for interesting discussions during the development of this work. The computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) and the Swedish National Infrastructure for Computing (SNIC) at NSC partially funded by the Swedish Research Council through grant agreements no. 2022-06725 and no. 2018-05973. This work is supported by the Swedish Research Council (project no. 2020-00452). SMC is supported by the US Department of Energy, Office of Science, Office of Nuclear Physics, Early Career Research Program under Award Number DE-SC0015904. Software. FLASH (Fryxell et al. 2000), NuLib (O’Connor 2015), Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), yt (Turk et al. 2011).

References

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016a, Phys. Rev. Lett., 116, 061102 [Google Scholar]
Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016b, Phys. Rev. Lett., 116, 241103 [Google Scholar]
Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016c, Phys. Rev. X, 6, 041015 [Google Scholar]
Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, ApJ, 851, L35 [Google Scholar]
Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017b, Phys. Rev. Lett., 119, 141101 [CrossRef] [Google Scholar]
Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017c, Phys. Rev. Lett., 119, 161101 [Google Scholar]
Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017d, Phys. Rev. Lett., 118, 221101 [CrossRef] [Google Scholar]
Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, Phys. Rev. X, 9, 031040 [Google Scholar]
Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2020a, ApJ, 892, L3 [NASA ADS] [CrossRef] [Google Scholar]
Abbott, R., Abbott, T. D., Abraham, S., et al. 2020b, Phys. Rev. Lett., 125, 101102 [Google Scholar]
Abbott, R., Abbott, T. D., Abraham, S., et al. 2020c, Phys. Rev. D, 102, 043015 [NASA ADS] [CrossRef] [Google Scholar]
Abbott, R., Abbott, T. D., Abraham, S., et al. 2020d, ApJ, 896, L44 [Google Scholar]
Abbott, R., Abbott, T. D., Abraham, S., et al. 2021a, ApJ, 915, L5 [NASA ADS] [CrossRef] [Google Scholar]
Abbott, R., Abbott, T. D., Abraham, S., et al. 2021b, Phys. Rev. X, 11, 021053 [Google Scholar]
Abbott, R., Abbott, T. D., Acernese, F., et al. 2023, Phys. Rev. X, 13, 041039 [Google Scholar]
Abbott, R., Abbott, T. D., Acernese, F., et al. 2024, Phys. Rev. D, 109, 022001 [NASA ADS] [CrossRef] [Google Scholar]
Abdikamalov, E., Burrows, A., Ott, C. D., et al. 2012, ApJ, 755, 111 [NASA ADS] [CrossRef] [Google Scholar]
Akaho, R., Harada, A., Nagakura, H., et al. 2021, ApJ, 909, 210 [NASA ADS] [CrossRef] [Google Scholar]
Akaho, R., Harada, A., Nagakura, H., et al. 2023, ApJ, 944, 60 [NASA ADS] [CrossRef] [Google Scholar]
Akaho, R., Liu, J., Nagakura, H., Zaizen, M., & Yamada, S. 2024, Phys. Rev. D, 109, 023012 [NASA ADS] [CrossRef] [Google Scholar]
Andresen, H. 2017, PhD Thesis, Munich, Tech. U., Germany [Google Scholar]
Andresen, H., Müller, B., Müller, E., & Janka, H.-T. 2017, MNRAS, 468, 2032 [NASA ADS] [CrossRef] [Google Scholar]
Andresen, H., Müller, E., Janka, H. T., et al. 2019, MNRAS, 486, 2238 [NASA ADS] [CrossRef] [Google Scholar]
Ardevol-Pulpillo, R., Janka, H. T., Just, O., & Bauswein, A. 2019, MNRAS, 485, 4754 [CrossRef] [Google Scholar]
Arnett, W. D., & Meakin, C. 2011, ApJ, 733, 78 [NASA ADS] [CrossRef] [Google Scholar]
Bacchini, C., Fraternali, F., Pezzulli, G., et al. 2023, in Resolving the Rise and Fall of Star Formation in Galaxies, 373, eds. T. Wong, & W.-T. Kim, 199 [Google Scholar]
Bhattacharyya, S., & Dasgupta, B. 2021, Phys. Rev. Lett., 126, 061302 [NASA ADS] [CrossRef] [Google Scholar]
Blanchet, L., Damour, T., & Schäfer, G. 1990, MNRAS, 242, 289 [NASA ADS] [CrossRef] [Google Scholar]
Bruenn, S. W., Blondin, J. M., Hix, W. R., et al. 2020, ApJS, 248, 11 [CrossRef] [Google Scholar]
Buellet, A. C., Foglizzo, T., Guilet, J., & Abdikamalov, E. 2023, Astron. Astrophys., 674, A205 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bugli, M., Guilet, J., & Obergaulinger, M. 2021, MNRAS, 507, 443 [NASA ADS] [CrossRef] [Google Scholar]
Bugli, M., Guilet, J., Foglizzo, T., & Obergaulinger, M. 2023, MNRAS, 520, 5622 [NASA ADS] [CrossRef] [Google Scholar]
Burrows, A., Young, T., Young, T., et al. 2000, ApJ, 539, 865 [NASA ADS] [CrossRef] [Google Scholar]
Burrows, A., Reddy, S., & Thompson, T. A. 2006, Nucl. Phys. A, 777, 356 [Google Scholar]
Burrows, A., Radice, D., & Vartanyan, D. 2019, MNRAS, 485, 3153 [NASA ADS] [CrossRef] [Google Scholar]
Burrows, A., Radice, D., Vartanyan, D., et al. 2020, MNRAS, 491, 2715 [NASA ADS] [CrossRef] [Google Scholar]
Capozzi, F., Raffelt, G., & Stirner, T. 2019, JCAP, 09, 002 [CrossRef] [Google Scholar]
Capozzi, F., Chakraborty, M., Chakraborty, S., & Sen, M. 2020, Phys. Rev. Lett., 125, 251801 [NASA ADS] [CrossRef] [Google Scholar]
Capozzi, F., Chakraborty, M., Chakraborty, S., & Sen, M. 2022, Phys. Rev. D, 106, 083011 [NASA ADS] [CrossRef] [Google Scholar]
Cardall, C. Y., Endeve, E., & Mezzacappa, A. 2013, Phys. Rev. D, 87, 103004 [NASA ADS] [CrossRef] [Google Scholar]
Chakraborty, M., & Chakraborty, S. 2020, JCAP, 01, 005 [CrossRef] [Google Scholar]
Chakraborty, S., Hansen, R. S., Izaguirre, I., & Raffelt, G. 2016, JCAP, 03, 042 [CrossRef] [Google Scholar]
Cornelius, M., Shalgar, S., & Tamborra, I. 2024, JCAP, 02, 038 [CrossRef] [Google Scholar]
Couch, S. M. 2013, ApJ, 765, 29 [CrossRef] [Google Scholar]
Couch, S. M., & O’Connor, E. P. 2014, ApJ, 785, 123 [NASA ADS] [CrossRef] [Google Scholar]
Couch, S. M., & Ott, C. D. 2015, ApJ, 799, 5 [CrossRef] [Google Scholar]
Cristini, A., Meakin, C., Hirschi, R., et al. 2017, MNRAS, 471, 279 [NASA ADS] [CrossRef] [Google Scholar]
Curtis, S., Bosch, P., Mösta, P., et al. 2024, ApJ, 961, L26 [NASA ADS] [CrossRef] [Google Scholar]
da Silva Schneider, A., O’Connor, E., Granqvist, E., Betranhandy, A., & Couch, S. M. 2020, ApJ, 894, 4 [NASA ADS] [CrossRef] [Google Scholar]
Dedin Neto, P., Tamborra, I., & Shalgar, S. 2023, arXiv e-prints [arXiv:2312.06556] [Google Scholar]
Edmunds, M. G. 2017, Supernovae and the Chemical Evolution of Galaxies, eds. A. W. Alsabti & P. Murdin (Cham: Springer International Publishing), 2455 [Google Scholar]
Eggenberger Andersen, O., Zha, S., da Silva Schneider, A., et al. 2021, ApJ, 923, 201 [NASA ADS] [CrossRef] [Google Scholar]
Ehring, J., Abbar, S., Janka, H.-T., Raffelt, G., & Tamborra, I. 2023, Phys. Rev. Lett., 131, 061401 [NASA ADS] [CrossRef] [Google Scholar]
Finn, L. S. 1989, in Frontiers in Numerical Relativity, eds. C. R. Evans, L. S. Finn, & D. W. Hobill (Cambridge (UK): Cambridge University Press), 126 [Google Scholar]
Fischer, T., Hempel, M., Sagert, I., Suwa, Y., & Schaffner-Bielich, J. 2014, Eur. Phys. J. A, 50, 46 [NASA ADS] [CrossRef] [Google Scholar]
Foucart, F. 2023, Living Rev. Computat. Astrophys., 9, 1 [NASA ADS] [CrossRef] [Google Scholar]
Foucart, F., O’Connor, E., Roberts, L., et al. 2015, Phys. Rev. D, 91, 124021 [CrossRef] [Google Scholar]
Foucart, F., O’Connor, E., Roberts, L., et al. 2016, Phys. Rev. D, 94, 123016 [NASA ADS] [CrossRef] [Google Scholar]
Fryer, C. L., & Warren, M. S. 2004, ApJ, 601, 391 [NASA ADS] [CrossRef] [Google Scholar]
Fryxell, B., Olson, K., Ricker, P., et al. 2000, The ApJS, 131, 273 [NASA ADS] [CrossRef] [Google Scholar]
Fujibayashi, S., Wanajo, S., Kiuchi, K., et al. 2020, ApJ, 901, 122 [NASA ADS] [CrossRef] [Google Scholar]
Fujibayashi, S., Kiuchi, K., Wanajo, S., et al. 2023, ApJ, 942, 39 [NASA ADS] [CrossRef] [Google Scholar]
Fukushima, T. 2015, Appl. Math. Comput., 259, 708 [Google Scholar]
Harada, A., Nagakura, H., Iwakami, W., et al. 2020, ApJ, 902, 150 [CrossRef] [Google Scholar]
Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
Hayashi, K., Fujibayashi, S., Kiuchi, K., et al. 2022, Phys. Rev. D, 106, 023008 [NASA ADS] [CrossRef] [Google Scholar]
Hempel, M., Fischer, T., Schaffner-Bielich, J., & Liebendörfer, M. 2012, ApJ, 748, 70 [CrossRef] [Google Scholar]
Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
Iwakami, W., Okawa, H., Nagakura, H., et al. 2020, ApJ, 903, 82 [NASA ADS] [CrossRef] [Google Scholar]
Izaguirre, I., Raffelt, G., & Tamborra, I. 2017, Phys. Rev. Lett., 118, 021101 [NASA ADS] [CrossRef] [Google Scholar]
Jardine, R., Powell, J., & Müller, B. 2022, MNRAS, 510, 5535 [NASA ADS] [CrossRef] [Google Scholar]
Johns, L. 2023, Phys. Rev. Lett., 130, 191001 [NASA ADS] [CrossRef] [Google Scholar]
Johns, L., Nagakura, H., Fuller, G. M., & Burrows, A. 2020, Phys. Rev. D, 101, 043009 [NASA ADS] [CrossRef] [Google Scholar]
Jones, S., Andrassy, R., Sandalski, S., et al. 2017, MNRAS, 465, 2991 [NASA ADS] [CrossRef] [Google Scholar]
Just, O., Obergaulinger, M., & Janka, H. T. 2015, MNRAS, 453, 3386 [Google Scholar]
Just, O., Bollig, R., Janka, H. T., et al. 2018, MNRAS, 481, 4786 [NASA ADS] [CrossRef] [Google Scholar]
Kato, C., Nagakura, H., Hori, Y., & Yamada, S. 2020, ApJ, 897, 43 [CrossRef] [Google Scholar]
Kiuchi, K., Fujibayashi, S., Hayashi, K., et al. 2023, Phys. Rev. Lett., 131, 011401 [NASA ADS] [CrossRef] [Google Scholar]
Kiuchi, K., Reboul-Salze, A., Shibata, M., & Sekiguchi, Y. 2024, Nat. Astron., 8, 298 [Google Scholar]
Kotake, K., Iwakami-Nakano, W., & Ohnishi, N. 2011, ApJ, 736, 124 [NASA ADS] [CrossRef] [Google Scholar]
Kotake, K., Takiwaki, T., Fischer, T., Nakamura, K., & Martínez-Pinedo, G. 2018, ApJ, 853, 170 [CrossRef] [Google Scholar]
Kuroda, T., Kotake, K., Takiwaki, T., & Thielemann, F.-K. 2018, MNRAS, 477, L80 [NASA ADS] [CrossRef] [Google Scholar]
Kuroda, T., Fischer, T., Takiwaki, T., & Kotake, K. 2022, ApJ, 924, 38 [NASA ADS] [CrossRef] [Google Scholar]
Liebendörfer, M., Messer, O. E. B., Mezzacappa, A., et al. 2004, ApJS, 150, 263 [CrossRef] [Google Scholar]
Lindquist, R. W. 1966, Ann. Phys., 37, 487 [NASA ADS] [CrossRef] [Google Scholar]
Liu, J., Akaho, R., Ito, A., et al. 2023, Phys. Rev. D, 108, 123024 [CrossRef] [Google Scholar]
Marek, A., Dimmelmeier, H., Janka, H. T., Müller, E., & Buras, R. 2006, A&A, 445, 273 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Martin, J. D., Carlson, J., Cirigliano, V., & Duan, H. 2021, Phys. Rev. D, 103, 063001 [CrossRef] [Google Scholar]
Matsumoto, J., Asahina, Y., Takiwaki, T., Kotake, K., & Takahashi, H. R. 2022, MNRAS, 516, 1752 [NASA ADS] [CrossRef] [Google Scholar]
Melson, T., Janka, H.-T., Bollig, R., et al. 2015a, ApJ, 808, L42 [NASA ADS] [CrossRef] [Google Scholar]
Melson, T., Janka, H.-T., & Marek, A. 2015b, ApJ, 801, L24 [NASA ADS] [CrossRef] [Google Scholar]
Mezzacappa, A., Endeve, E., Bronson Messer, O. E., & Bruenn, S. W. 2020, Liv. Rev. Comput. Astrophys., 6, 4 [CrossRef] [Google Scholar]
Mezzacappa, A., Marronetti, P., Landfield, R. E., et al. 2023, Phys. Rev. D, 107, 043008 [NASA ADS] [CrossRef] [Google Scholar]
Mori, M., Suwa, Y., & Takiwaki, T. 2023, Phys. Rev. D, 107, 083015 [NASA ADS] [CrossRef] [Google Scholar]
Morozova, V., Radice, D., Burrows, A., & Vartanyan, D. 2018, ApJ, 861, 10 [CrossRef] [Google Scholar]
Mueller, B., Janka, H.-T., & Marek, A. 2013, ApJ, 766, 43 [CrossRef] [Google Scholar]
Müller, B., & Janka, H. T. 2015, MNRAS, 448, 2141 [CrossRef] [Google Scholar]
Müller, B., Janka, H.-T., & Marek, A. 2013, ApJ, 766, 43 [CrossRef] [Google Scholar]
Müller, B., Viallet, M., Heger, A., & Janka, H.-T. 2016, ApJ, 833, 124 [CrossRef] [Google Scholar]
Murphy, J. W., Ott, C. D., & Burrows, A. 2009, ApJ, 707, 1173 [NASA ADS] [CrossRef] [Google Scholar]
Nagakura, H., Sumiyoshi, K., & Yamada, S. 2014, ApJS, 214, 16 [CrossRef] [Google Scholar]
Nagakura, H., Iwakami, W., Furusawa, S., et al. 2018, ApJ, 854, 136 [CrossRef] [Google Scholar]
Nagakura, H., Sumiyoshi, K., & Yamada, S. 2019, ApJ, 878, 160 [NASA ADS] [CrossRef] [Google Scholar]
Obergaulinger, M., & Aloy, M.-A. 2021, MNRAS, 503, 4942 [NASA ADS] [CrossRef] [Google Scholar]
O’Connor, E. 2015, ApJS, 219, 24 [CrossRef] [Google Scholar]
O’Connor, E. P., & Couch, S. M. 2018a, ApJ, 865, 81 [CrossRef] [Google Scholar]
O’Connor, E. P., & Couch, S. M. 2018b, ApJ, 854, 63 [CrossRef] [Google Scholar]
O’Connor, E., Bollig, R., Burrows, A., et al. 2018, J. Phys. G Nucl. Phys., 45, 104001 [CrossRef] [Google Scholar]
Oertel, M., Hempel, M., Klähn, T., & Typel, S. 2017, Rev. Mod. Phys., 89, 015007 [Google Scholar]
Oohara, K.-i., Nakamura, T., & Shibata, M. 1997, Prog. Theor. Phys. Supp., 128, 183 [CrossRef] [Google Scholar]
Pan, K.-C., Liebendörfer, M., Hempel, M., & Thielemann, F.-K. 2016, ApJ, 817, 72 [NASA ADS] [CrossRef] [Google Scholar]
Pascal, A., Novak, J., & Oertel, M. 2022, MNRAS, 511, 356 [NASA ADS] [CrossRef] [Google Scholar]
Radice, D., & Bernuzzi, S. 2023, ApJ, 959, 46 [NASA ADS] [CrossRef] [Google Scholar]
Radice, D., Morozova, V., Burrows, A., Vartanyan, D., & Nagakura, H. 2019, ApJ, 876, L9 [NASA ADS] [CrossRef] [Google Scholar]
Radice, D., Bernuzzi, S., Perego, A., & Haas, R. 2022, MNRAS, 512, 1499 [NASA ADS] [CrossRef] [Google Scholar]
Rahman, N., Just, O., & Janka, H. T. 2019, MNRAS, 490, 3545 [NASA ADS] [CrossRef] [Google Scholar]
Rahman, N., Janka, H.-T., Stockinger, G., & Woosley, S. 2022, MNRAS, 512, 4503 [NASA ADS] [CrossRef] [Google Scholar]
Reichert, M., Obergaulinger, M., Aloy, M. Á., et al. 2023, MNRAS, 518, 1557 [Google Scholar]
Richers, S., & Sen, M. 2022, Fast Flavor Transformations, eds. I. Tanihata, H. Toki, & T. Kajino (Springer), 1 [Google Scholar]
Richers, S., Nagakura, H., Ott, C. D., et al. 2017, ApJ, 847, 133 [NASA ADS] [CrossRef] [Google Scholar]
Roberts, L. F., Ott, C. D., Haas, R., et al. 2016, ApJ, 831, 98 [CrossRef] [Google Scholar]
Rodriguez, M. C., Ranea-Sandoval, I. F., Chirenti, C., & Radice, D. 2023, MNRAS, 523, 2236 [NASA ADS] [CrossRef] [Google Scholar]
Rosswog, S., & Liebendoerfer, M. 2003, MNRAS, 342, 673 [NASA ADS] [CrossRef] [Google Scholar]
Scheck, L., Kifonidis, K., Janka, H. T., & Mueller, E. 2006, Astron. Astrophys., 457, 963 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Schianchi, F., Gieg, H., Nedora, V., et al. 2024, Phys. Rev. D, 109, 044012 [CrossRef] [Google Scholar]
Schneider, A. S., Roberts, L. F., & Ott, C. D. 2017, Phys. Rev. C, 96, 065802 [Google Scholar]
Schneider, A. S., Roberts, L. F., Ott, C. D., & O’Connor, E. 2019, Phys. Rev. C, 100, 055802 [CrossRef] [Google Scholar]
Shalgar, S., & Tamborra, I. 2024, Phys. Rev. D, 109, 103011 [NASA ADS] [CrossRef] [Google Scholar]
Shibata, M., Kiuchi, K., Sekiguchi, Y.-i., & Suwa, Y. 2011, Prog. Theor. Phys., 125, 1255 [CrossRef] [Google Scholar]
Smith, M. C., Sijacki, D., & Shen, S. 2018, MNRAS, 478, 302 [NASA ADS] [CrossRef] [Google Scholar]
Sotani, H., & Takiwaki, T. 2016, Phys. Rev. D, 94, 044043 [NASA ADS] [CrossRef] [Google Scholar]
Sotani, H., Kuroda, T., Takiwaki, T., & Kotake, K. 2017, Phys. Rev. D, 96, 063005 [NASA ADS] [CrossRef] [Google Scholar]
Sotani, H., Kuroda, T., Takiwaki, T., & Kotake, K. 2019, Phys. Rev. D, 99, 123024 [CrossRef] [Google Scholar]
Sotani, H., Takiwaki, T., & Togashi, H. 2021, Phys. Rev. D, 104, 123009 [NASA ADS] [CrossRef] [Google Scholar]
Steiner, A. W., Hempel, M., & Fischer, T. 2013, ApJ, 774, 17 [NASA ADS] [CrossRef] [Google Scholar]
Sukhbold, T., Woosley, S. E., & Heger, A. 2018, ApJ, 860, 93 [NASA ADS] [CrossRef] [Google Scholar]
Suleiman, L., Oertel, M., & Mancini, M. 2023, Phys. Rev. C, 108, 035803 [NASA ADS] [CrossRef] [Google Scholar]
Sullivan, C., O’Connor, E., Zegers, R. G. T., Grubb, T., & Austin, S. M. 2016, ApJ, 816, 44 [NASA ADS] [CrossRef] [Google Scholar]
Sumiyoshi, K., Yamada, S., Suzuki, H., et al. 2005, ApJ, 629, 922 [NASA ADS] [CrossRef] [Google Scholar]
Summa, A., Janka, H. T., Melson, T., & Marek, A. 2018, ApJ, 852, 28 [NASA ADS] [CrossRef] [Google Scholar]
Takiwaki, T., & Kotake, K. 2018, MNRAS, 475, L91 [NASA ADS] [CrossRef] [Google Scholar]
Takiwaki, T., Kotake, K., & Suwa, Y. 2014, ApJ, 786, 83 [NASA ADS] [CrossRef] [Google Scholar]
Thielemann, F.-K., Isern, J., Perego, A., & von Ballmoos, P. 2018, Space Sci. Rev., 214, 62 [CrossRef] [Google Scholar]
Thorne, K. S. 1980, Rev. Mod. Phys., 52, 299 [NASA ADS] [CrossRef] [Google Scholar]
Torres-Forné, A., Cerdá-Durán, P., Passamonti, A., & Font, J. A. 2018, MNRAS, 474, 5272 [CrossRef] [Google Scholar]
Torres-Forné, A., Cerdá-Durán, P., Obergaulinger, M., Müller, B., & Font, J. A. 2019a, Phys. Rev. Lett., 123, 051102; Erratum: 2021, 127, 239901 [CrossRef] [Google Scholar]
Torres-Forné, A., Cerdá-Durán, P., Passamonti, A., Obergaulinger, M., & Font, J. A. 2019b, MNRAS, 482, 3967 [CrossRef] [Google Scholar]
Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9 [NASA ADS] [CrossRef] [Google Scholar]
Vartanyan, D., Burrows, A., Radice, D., Skinner, M. A., & Dolence, J. 2018, MNRAS, 477, 3091 [NASA ADS] [CrossRef] [Google Scholar]
Vartanyan, D., Burrows, A., Radice, D., Skinner, A. M., & Dolence, J. 2019, MNRAS, 482, 351 [NASA ADS] [CrossRef] [Google Scholar]
Vartanyan, D., Coleman, M. S. B., & Burrows, A. 2022, MNRAS, 510, 4689 [NASA ADS] [CrossRef] [Google Scholar]
Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
Wang, T., & Burrows, A. 2024, ApJ, 962, 71 [NASA ADS] [CrossRef] [Google Scholar]
Westernacher-Schneider, J. R., O’Connor, E., O’Sullivan, E., et al. 2019, Phys. Rev. D, 100, 123009 [NASA ADS] [CrossRef] [Google Scholar]
Wolfe, N. E., Fröhlich, C., Miller, J. M., Torres-Forné, A., & Cerdá-Durán, P. 2023, ApJ, 954, 161 [NASA ADS] [CrossRef] [Google Scholar]
Woosley, S. E., & Heger, A. 2007, Phys. Rep., 442, 269 [Google Scholar]
Xiong, Z., Wu, M.-R., Abbar, S., et al. 2023a, Phys. Rev. D, 108, 063003 [NASA ADS] [CrossRef] [Google Scholar]
Xiong, Z., Wu, M.-R., Martínez-Pinedo, G., et al. 2023b, Phys. Rev. D, 107, 083016 [CrossRef] [Google Scholar]
Yasin, H., Schäfer, S., Arcones, A., & Schwenk, A. 2020, Phys. Rev. Lett., 124, 092701 [NASA ADS] [CrossRef] [Google Scholar]
Zappa, F., Bernuzzi, S., Radice, D., & Perego, A. 2023, MNRAS, 520, 1481 [Google Scholar]

¹

These integrals can be numerically challenging and we are required to evaluate a large number of them for each time step. We followed the accurate and fast approach of Fukushima (2015) to compute the integrals.

²

https://github.com/evanoconnor/NuLib

³

The spike seen in the red curve in the middle panel is due to numerical issues arising in the integral in Eq. (30) when H^r is close to zero.

All Tables

Table 1

Neutrino luminosities and average neutrino energy in hydrostatic simulations of the post bounce phase of a 30 solar mass progenitor.

In the text

Table 2

Overview of the model names for our second set of 2D simulations.

In the text

All Figures

Fig. 1

Neutrino luminosities for the gray simulation compared to the simulations presented in O’Connor et al. (2018; semitransparent colored lines). The purple lines indicate the results from the simulation based on the original gray transport from Foucart et al. (2016), and the black lines represent the result from the simulation using the improvements detailed in this work. The luminosity of the electron neutrinos (solid lines) and electron anti-neutrinos (dot-dashed lines) are shown in the left panel, and the luminosity of the heavy-lepton neutrinos is shown in the right panel (dashed lines). Time is given in seconds after bounce. The inset in the left panel shows a zoom-in of the peak centered around ~0.l s after bounce.

In the text

Fig. 2

Average neutrino energies for the gray simulation (bright lines) compared to the simulations presented in O’Connor et al. (2018; semitransparent colored lines). The purple lines indicate the results from the simulation based on the original gray transport from Foucart et al. (2016), and the black lines represent the result from the simulation using the improvements detailed in this work. The average energy of the electron neutrinos (solid lines) and electron anti-neutrinos (dot-dashed lines) are shown in the left panel, and the average energies of the heavy-lepton neutrinos are shown in the right panel (dashed lines). Time is given in seconds after bounce.

In the text

Fig. 3

Shock radius (solid lines) and PNS radius (dashed lines) for the gray simulation (purple lines) compared to the simulations presented in O’Connor et al. (2018; semitransparent colored lines). The purple lines indicate the results from the simulation based on the original gray transport from Foucart et al. (2016) and the black lines represent the result from the simulation using the improvements detailed in this work. Time is given in seconds after bounce. The PNS radius is defined to be where the density drops below 10¹¹ g cm⁻³.

In the text

Fig. 4

Energy-averaged absorption opacity as a function of radius for a hydrostatic evolution of a 30 solar mass progenitor using the SFHo EOS. The dashed red lines show the opacities calculated with the neutrino spectrum from the energy-dependent neutrino transport (Eq. (29)). The dashed green lines show the opacities calculated according to Eq. (19). The solid blue lines represent the value according to Eq. (32). The panels show results for electron neutrinos (top), electron anti-neutrinos (middle), and heavy-lepton neutrinos (bottom). The opacities for the heavy-lepton neutrinos refer to effective opacities, which incorporate pair-production processes in a simplified way.

In the text

Fig. 5

Sum of the energy-averaged absorption and energy-averaged scattering opacities as a function of radius for a hydrostatic simulation of a 30 solar mass progenitor using the SFHo EOS. The dashed red lines show the opacities calculated with the neutrino spectrum from the energy-dependent neutrino transport (Eq. (30)). The dashed green lines show the opacities calculated according to Eq. (19). The blue lines represent the value according to Eq. (33). The panels show results for electron neutrinos (top), electron anti-neutrinos (middle), and heavy-lepton neutrinos (bottom). The opacities for the heavy-lepton neutrinos refer to effective opacities, which incorporate pair-production processes in a simplified way.

In the text

Fig. 6

Average shock radius (left column) and the PNS radius (right column) for simulations based on five progenitors with ZAMS masses of 11, 12, 13, 14, and 15 solar masses. We performed three simulations for each progenitor, using different neutrino transport methods. The dashed red line corresponds to simulations with fully energy-dependent transport, the green lines indicates simulations using the gray scheme of Foucart et al. (2016), and the blue lines show results from simulations using the gray transport proposed in this work. Each row corresponds to a given progenitor; the progenitor model is indicated in the top left corner of each row. The simulations were performed with the SFHo EOS.

In the text

Fig. 7

Average shock radius (left column) and the PNS radius (right column) for simulations based on one progenitor with a ZAMS mass of 20 solar masses. For each EOS we performed three simulations, one for each neutrino transport method. The dashed red line corresponds to simulations with fully energy-dependent transport, the green lines indicate simulations using the gray scheme of Foucart et al. (2016), and the blue lines show results from simulations using the gray transport proposed in this work. Each row corresponds to a given version of the SRO EOS; EOS parameters (see text) are indicated by the model name in the top left corner of each row.

In the text

Fig. 8

Heating rate as a function of time for models s11, s12, s13, s14, and s15. The blue curves indicate the results from simulations with our gray scheme, the green curves show the results from simulations using the neutrino transport of Foucart et al. (201 ), and the red curves correspond to the results from simulations with energy-dependent neutrino transport. Each row correspond to one model; the model names are shown in the bottom left corner of each panel.

In the text

Fig. 9

GWs from two simulations of the s12 progenitor. The right column shows the signals from the simulation with energy-dependent neutrino transport and the left column shows the corresponding simulation performed with our gray neutrino transport. The top row shows the strain as a function of time and the bottom row shows the square of the short time Fourier transform calculated with scipy.signal.stft (Virtanen et al. 2020). The color scale is logarithmic, and the plots were normalized by a common factor.

In the text

Fig. 10

Neutrino luminosities (top panel) and average energies (bottom panel) from two simulations of the s12 progenitor. The red lines indicate the results from the simulation with energy-dependent neutrino transport. The blue lines show the results from the simulation performed with our gray neutrino transport. The solid lines indicate electron neutrinos, dashed lines show electron anti-neutrinos, and dot-dashed lines correspond to heavy-lepton neutrinos.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.

[1] Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016a, Phys. Rev. Lett., 116, 061102 [Google Scholar]

[2] Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016b, Phys. Rev. Lett., 116, 241103 [Google Scholar]

[3] Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016c, Phys. Rev. X, 6, 041015 [Google Scholar]

[4] Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, ApJ, 851, L35 [Google Scholar]

[5] Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017b, Phys. Rev. Lett., 119, 141101 [CrossRef] [Google Scholar]

[6] Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017c, Phys. Rev. Lett., 119, 161101 [Google Scholar]

[7] Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017d, Phys. Rev. Lett., 118, 221101 [CrossRef] [Google Scholar]

[8] Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, Phys. Rev. X, 9, 031040 [Google Scholar]

[9] Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2020a, ApJ, 892, L3 [NASA ADS] [CrossRef] [Google Scholar]

[10] Abbott, R., Abbott, T. D., Abraham, S., et al. 2020b, Phys. Rev. Lett., 125, 101102 [Google Scholar]

[11] Abbott, R., Abbott, T. D., Abraham, S., et al. 2020c, Phys. Rev. D, 102, 043015 [NASA ADS] [CrossRef] [Google Scholar]

[12] Abbott, R., Abbott, T. D., Abraham, S., et al. 2020d, ApJ, 896, L44 [Google Scholar]

[13] Abbott, R., Abbott, T. D., Abraham, S., et al. 2021a, ApJ, 915, L5 [NASA ADS] [CrossRef] [Google Scholar]

[14] Abbott, R., Abbott, T. D., Abraham, S., et al. 2021b, Phys. Rev. X, 11, 021053 [Google Scholar]

[15] Abbott, R., Abbott, T. D., Acernese, F., et al. 2023, Phys. Rev. X, 13, 041039 [Google Scholar]

[16] Abbott, R., Abbott, T. D., Acernese, F., et al. 2024, Phys. Rev. D, 109, 022001 [NASA ADS] [CrossRef] [Google Scholar]

[17] Abdikamalov, E., Burrows, A., Ott, C. D., et al. 2012, ApJ, 755, 111 [NASA ADS] [CrossRef] [Google Scholar]

[18] Akaho, R., Harada, A., Nagakura, H., et al. 2021, ApJ, 909, 210 [NASA ADS] [CrossRef] [Google Scholar]

[19] Akaho, R., Harada, A., Nagakura, H., et al. 2023, ApJ, 944, 60 [NASA ADS] [CrossRef] [Google Scholar]

[20] Akaho, R., Liu, J., Nagakura, H., Zaizen, M., & Yamada, S. 2024, Phys. Rev. D, 109, 023012 [NASA ADS] [CrossRef] [Google Scholar]

[21] Andresen, H. 2017, PhD Thesis, Munich, Tech. U., Germany [Google Scholar]

[22] Andresen, H., Müller, B., Müller, E., & Janka, H.-T. 2017, MNRAS, 468, 2032 [NASA ADS] [CrossRef] [Google Scholar]

[23] Andresen, H., Müller, E., Janka, H. T., et al. 2019, MNRAS, 486, 2238 [NASA ADS] [CrossRef] [Google Scholar]

[24] Ardevol-Pulpillo, R., Janka, H. T., Just, O., & Bauswein, A. 2019, MNRAS, 485, 4754 [CrossRef] [Google Scholar]

[25] Arnett, W. D., & Meakin, C. 2011, ApJ, 733, 78 [NASA ADS] [CrossRef] [Google Scholar]

[26] Bacchini, C., Fraternali, F., Pezzulli, G., et al. 2023, in Resolving the Rise and Fall of Star Formation in Galaxies, 373, eds. T. Wong, & W.-T. Kim, 199 [Google Scholar]

[27] Bhattacharyya, S., & Dasgupta, B. 2021, Phys. Rev. Lett., 126, 061302 [NASA ADS] [CrossRef] [Google Scholar]

[28] Blanchet, L., Damour, T., & Schäfer, G. 1990, MNRAS, 242, 289 [NASA ADS] [CrossRef] [Google Scholar]

[29] Bruenn, S. W., Blondin, J. M., Hix, W. R., et al. 2020, ApJS, 248, 11 [CrossRef] [Google Scholar]

[30] Buellet, A. C., Foglizzo, T., Guilet, J., & Abdikamalov, E. 2023, Astron. Astrophys., 674, A205 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[31] Bugli, M., Guilet, J., & Obergaulinger, M. 2021, MNRAS, 507, 443 [NASA ADS] [CrossRef] [Google Scholar]

[32] Bugli, M., Guilet, J., Foglizzo, T., & Obergaulinger, M. 2023, MNRAS, 520, 5622 [NASA ADS] [CrossRef] [Google Scholar]

[33] Burrows, A., Young, T., Young, T., et al. 2000, ApJ, 539, 865 [NASA ADS] [CrossRef] [Google Scholar]

[34] Burrows, A., Reddy, S., & Thompson, T. A. 2006, Nucl. Phys. A, 777, 356 [Google Scholar]

[35] Burrows, A., Radice, D., & Vartanyan, D. 2019, MNRAS, 485, 3153 [NASA ADS] [CrossRef] [Google Scholar]

[36] Burrows, A., Radice, D., Vartanyan, D., et al. 2020, MNRAS, 491, 2715 [NASA ADS] [CrossRef] [Google Scholar]

[37] Capozzi, F., Raffelt, G., & Stirner, T. 2019, JCAP, 09, 002 [CrossRef] [Google Scholar]

[38] Capozzi, F., Chakraborty, M., Chakraborty, S., & Sen, M. 2020, Phys. Rev. Lett., 125, 251801 [NASA ADS] [CrossRef] [Google Scholar]

[39] Capozzi, F., Chakraborty, M., Chakraborty, S., & Sen, M. 2022, Phys. Rev. D, 106, 083011 [NASA ADS] [CrossRef] [Google Scholar]

[40] Cardall, C. Y., Endeve, E., & Mezzacappa, A. 2013, Phys. Rev. D, 87, 103004 [NASA ADS] [CrossRef] [Google Scholar]

[41] Chakraborty, M., & Chakraborty, S. 2020, JCAP, 01, 005 [CrossRef] [Google Scholar]

[42] Chakraborty, S., Hansen, R. S., Izaguirre, I., & Raffelt, G. 2016, JCAP, 03, 042 [CrossRef] [Google Scholar]

[43] Cornelius, M., Shalgar, S., & Tamborra, I. 2024, JCAP, 02, 038 [CrossRef] [Google Scholar]

[44] Couch, S. M. 2013, ApJ, 765, 29 [CrossRef] [Google Scholar]

[45] Couch, S. M., & O’Connor, E. P. 2014, ApJ, 785, 123 [NASA ADS] [CrossRef] [Google Scholar]

[46] Couch, S. M., & Ott, C. D. 2015, ApJ, 799, 5 [CrossRef] [Google Scholar]

[47] Cristini, A., Meakin, C., Hirschi, R., et al. 2017, MNRAS, 471, 279 [NASA ADS] [CrossRef] [Google Scholar]

[48] Curtis, S., Bosch, P., Mösta, P., et al. 2024, ApJ, 961, L26 [NASA ADS] [CrossRef] [Google Scholar]

[49] da Silva Schneider, A., O’Connor, E., Granqvist, E., Betranhandy, A., & Couch, S. M. 2020, ApJ, 894, 4 [NASA ADS] [CrossRef] [Google Scholar]

[50] Dedin Neto, P., Tamborra, I., & Shalgar, S. 2023, arXiv e-prints [arXiv:2312.06556] [Google Scholar]

[51] Edmunds, M. G. 2017, Supernovae and the Chemical Evolution of Galaxies, eds. A. W. Alsabti & P. Murdin (Cham: Springer International Publishing), 2455 [Google Scholar]

[52] Eggenberger Andersen, O., Zha, S., da Silva Schneider, A., et al. 2021, ApJ, 923, 201 [NASA ADS] [CrossRef] [Google Scholar]

[53] Ehring, J., Abbar, S., Janka, H.-T., Raffelt, G., & Tamborra, I. 2023, Phys. Rev. Lett., 131, 061401 [NASA ADS] [CrossRef] [Google Scholar]

[54] Finn, L. S. 1989, in Frontiers in Numerical Relativity, eds. C. R. Evans, L. S. Finn, & D. W. Hobill (Cambridge (UK): Cambridge University Press), 126 [Google Scholar]

[55] Fischer, T., Hempel, M., Sagert, I., Suwa, Y., & Schaffner-Bielich, J. 2014, Eur. Phys. J. A, 50, 46 [NASA ADS] [CrossRef] [Google Scholar]

[56] Foucart, F. 2023, Living Rev. Computat. Astrophys., 9, 1 [NASA ADS] [CrossRef] [Google Scholar]

[57] Foucart, F., O’Connor, E., Roberts, L., et al. 2015, Phys. Rev. D, 91, 124021 [CrossRef] [Google Scholar]

[58] Foucart, F., O’Connor, E., Roberts, L., et al. 2016, Phys. Rev. D, 94, 123016 [NASA ADS] [CrossRef] [Google Scholar]

[59] Fryer, C. L., & Warren, M. S. 2004, ApJ, 601, 391 [NASA ADS] [CrossRef] [Google Scholar]

[60] Fryxell, B., Olson, K., Ricker, P., et al. 2000, The ApJS, 131, 273 [NASA ADS] [CrossRef] [Google Scholar]

[61] Fujibayashi, S., Wanajo, S., Kiuchi, K., et al. 2020, ApJ, 901, 122 [NASA ADS] [CrossRef] [Google Scholar]

[62] Fujibayashi, S., Kiuchi, K., Wanajo, S., et al. 2023, ApJ, 942, 39 [NASA ADS] [CrossRef] [Google Scholar]

[63] Fukushima, T. 2015, Appl. Math. Comput., 259, 708 [Google Scholar]

[64] Harada, A., Nagakura, H., Iwakami, W., et al. 2020, ApJ, 902, 150 [CrossRef] [Google Scholar]

[65] Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]

[66] Hayashi, K., Fujibayashi, S., Kiuchi, K., et al. 2022, Phys. Rev. D, 106, 023008 [NASA ADS] [CrossRef] [Google Scholar]

[67] Hempel, M., Fischer, T., Schaffner-Bielich, J., & Liebendörfer, M. 2012, ApJ, 748, 70 [CrossRef] [Google Scholar]

[68] Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]

[69] Iwakami, W., Okawa, H., Nagakura, H., et al. 2020, ApJ, 903, 82 [NASA ADS] [CrossRef] [Google Scholar]

[70] Izaguirre, I., Raffelt, G., & Tamborra, I. 2017, Phys. Rev. Lett., 118, 021101 [NASA ADS] [CrossRef] [Google Scholar]

[71] Jardine, R., Powell, J., & Müller, B. 2022, MNRAS, 510, 5535 [NASA ADS] [CrossRef] [Google Scholar]

[72] Johns, L. 2023, Phys. Rev. Lett., 130, 191001 [NASA ADS] [CrossRef] [Google Scholar]

[73] Johns, L., Nagakura, H., Fuller, G. M., & Burrows, A. 2020, Phys. Rev. D, 101, 043009 [NASA ADS] [CrossRef] [Google Scholar]

[74] Jones, S., Andrassy, R., Sandalski, S., et al. 2017, MNRAS, 465, 2991 [NASA ADS] [CrossRef] [Google Scholar]

[75] Just, O., Obergaulinger, M., & Janka, H. T. 2015, MNRAS, 453, 3386 [Google Scholar]

[76] Just, O., Bollig, R., Janka, H. T., et al. 2018, MNRAS, 481, 4786 [NASA ADS] [CrossRef] [Google Scholar]

[77] Kato, C., Nagakura, H., Hori, Y., & Yamada, S. 2020, ApJ, 897, 43 [CrossRef] [Google Scholar]

[78] Kiuchi, K., Fujibayashi, S., Hayashi, K., et al. 2023, Phys. Rev. Lett., 131, 011401 [NASA ADS] [CrossRef] [Google Scholar]

[79] Kiuchi, K., Reboul-Salze, A., Shibata, M., & Sekiguchi, Y. 2024, Nat. Astron., 8, 298 [Google Scholar]

[80] Kotake, K., Iwakami-Nakano, W., & Ohnishi, N. 2011, ApJ, 736, 124 [NASA ADS] [CrossRef] [Google Scholar]

[81] Kotake, K., Takiwaki, T., Fischer, T., Nakamura, K., & Martínez-Pinedo, G. 2018, ApJ, 853, 170 [CrossRef] [Google Scholar]

[82] Kuroda, T., Kotake, K., Takiwaki, T., & Thielemann, F.-K. 2018, MNRAS, 477, L80 [NASA ADS] [CrossRef] [Google Scholar]

[83] Kuroda, T., Fischer, T., Takiwaki, T., & Kotake, K. 2022, ApJ, 924, 38 [NASA ADS] [CrossRef] [Google Scholar]

[84] Liebendörfer, M., Messer, O. E. B., Mezzacappa, A., et al. 2004, ApJS, 150, 263 [CrossRef] [Google Scholar]

[85] Lindquist, R. W. 1966, Ann. Phys., 37, 487 [NASA ADS] [CrossRef] [Google Scholar]

[86] Liu, J., Akaho, R., Ito, A., et al. 2023, Phys. Rev. D, 108, 123024 [CrossRef] [Google Scholar]

[87] Marek, A., Dimmelmeier, H., Janka, H. T., Müller, E., & Buras, R. 2006, A&A, 445, 273 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[88] Martin, J. D., Carlson, J., Cirigliano, V., & Duan, H. 2021, Phys. Rev. D, 103, 063001 [CrossRef] [Google Scholar]

[89] Matsumoto, J., Asahina, Y., Takiwaki, T., Kotake, K., & Takahashi, H. R. 2022, MNRAS, 516, 1752 [NASA ADS] [CrossRef] [Google Scholar]

[90] Melson, T., Janka, H.-T., Bollig, R., et al. 2015a, ApJ, 808, L42 [NASA ADS] [CrossRef] [Google Scholar]

[91] Melson, T., Janka, H.-T., & Marek, A. 2015b, ApJ, 801, L24 [NASA ADS] [CrossRef] [Google Scholar]

[92] Mezzacappa, A., Endeve, E., Bronson Messer, O. E., & Bruenn, S. W. 2020, Liv. Rev. Comput. Astrophys., 6, 4 [CrossRef] [Google Scholar]

[93] Mezzacappa, A., Marronetti, P., Landfield, R. E., et al. 2023, Phys. Rev. D, 107, 043008 [NASA ADS] [CrossRef] [Google Scholar]

[94] Mori, M., Suwa, Y., & Takiwaki, T. 2023, Phys. Rev. D, 107, 083015 [NASA ADS] [CrossRef] [Google Scholar]

[95] Morozova, V., Radice, D., Burrows, A., & Vartanyan, D. 2018, ApJ, 861, 10 [CrossRef] [Google Scholar]

[96] Mueller, B., Janka, H.-T., & Marek, A. 2013, ApJ, 766, 43 [CrossRef] [Google Scholar]

[97] Müller, B., & Janka, H. T. 2015, MNRAS, 448, 2141 [CrossRef] [Google Scholar]

[98] Müller, B., Janka, H.-T., & Marek, A. 2013, ApJ, 766, 43 [CrossRef] [Google Scholar]

[99] Müller, B., Viallet, M., Heger, A., & Janka, H.-T. 2016, ApJ, 833, 124 [CrossRef] [Google Scholar]

[100] Murphy, J. W., Ott, C. D., & Burrows, A. 2009, ApJ, 707, 1173 [NASA ADS] [CrossRef] [Google Scholar]

[101] Nagakura, H., Sumiyoshi, K., & Yamada, S. 2014, ApJS, 214, 16 [CrossRef] [Google Scholar]

[102] Nagakura, H., Iwakami, W., Furusawa, S., et al. 2018, ApJ, 854, 136 [CrossRef] [Google Scholar]

[103] Nagakura, H., Sumiyoshi, K., & Yamada, S. 2019, ApJ, 878, 160 [NASA ADS] [CrossRef] [Google Scholar]

[104] Obergaulinger, M., & Aloy, M.-A. 2021, MNRAS, 503, 4942 [NASA ADS] [CrossRef] [Google Scholar]

[105] O’Connor, E. 2015, ApJS, 219, 24 [CrossRef] [Google Scholar]

[106] O’Connor, E. P., & Couch, S. M. 2018a, ApJ, 865, 81 [CrossRef] [Google Scholar]

[107] O’Connor, E. P., & Couch, S. M. 2018b, ApJ, 854, 63 [CrossRef] [Google Scholar]

[108] O’Connor, E., Bollig, R., Burrows, A., et al. 2018, J. Phys. G Nucl. Phys., 45, 104001 [CrossRef] [Google Scholar]

[109] Oertel, M., Hempel, M., Klähn, T., & Typel, S. 2017, Rev. Mod. Phys., 89, 015007 [Google Scholar]

[110] Oohara, K.-i., Nakamura, T., & Shibata, M. 1997, Prog. Theor. Phys. Supp., 128, 183 [CrossRef] [Google Scholar]

[111] Pan, K.-C., Liebendörfer, M., Hempel, M., & Thielemann, F.-K. 2016, ApJ, 817, 72 [NASA ADS] [CrossRef] [Google Scholar]

[112] Pascal, A., Novak, J., & Oertel, M. 2022, MNRAS, 511, 356 [NASA ADS] [CrossRef] [Google Scholar]

[113] Radice, D., & Bernuzzi, S. 2023, ApJ, 959, 46 [NASA ADS] [CrossRef] [Google Scholar]

[114] Radice, D., Morozova, V., Burrows, A., Vartanyan, D., & Nagakura, H. 2019, ApJ, 876, L9 [NASA ADS] [CrossRef] [Google Scholar]

[115] Radice, D., Bernuzzi, S., Perego, A., & Haas, R. 2022, MNRAS, 512, 1499 [NASA ADS] [CrossRef] [Google Scholar]

[116] Rahman, N., Just, O., & Janka, H. T. 2019, MNRAS, 490, 3545 [NASA ADS] [CrossRef] [Google Scholar]

[117] Rahman, N., Janka, H.-T., Stockinger, G., & Woosley, S. 2022, MNRAS, 512, 4503 [NASA ADS] [CrossRef] [Google Scholar]

[118] Reichert, M., Obergaulinger, M., Aloy, M. Á., et al. 2023, MNRAS, 518, 1557 [Google Scholar]

[119] Richers, S., & Sen, M. 2022, Fast Flavor Transformations, eds. I. Tanihata, H. Toki, & T. Kajino (Springer), 1 [Google Scholar]

[120] Richers, S., Nagakura, H., Ott, C. D., et al. 2017, ApJ, 847, 133 [NASA ADS] [CrossRef] [Google Scholar]

[121] Roberts, L. F., Ott, C. D., Haas, R., et al. 2016, ApJ, 831, 98 [CrossRef] [Google Scholar]

[122] Rodriguez, M. C., Ranea-Sandoval, I. F., Chirenti, C., & Radice, D. 2023, MNRAS, 523, 2236 [NASA ADS] [CrossRef] [Google Scholar]

[123] Rosswog, S., & Liebendoerfer, M. 2003, MNRAS, 342, 673 [NASA ADS] [CrossRef] [Google Scholar]

[124] Scheck, L., Kifonidis, K., Janka, H. T., & Mueller, E. 2006, Astron. Astrophys., 457, 963 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[125] Schianchi, F., Gieg, H., Nedora, V., et al. 2024, Phys. Rev. D, 109, 044012 [CrossRef] [Google Scholar]

[126] Schneider, A. S., Roberts, L. F., & Ott, C. D. 2017, Phys. Rev. C, 96, 065802 [Google Scholar]

[127] Schneider, A. S., Roberts, L. F., Ott, C. D., & O’Connor, E. 2019, Phys. Rev. C, 100, 055802 [CrossRef] [Google Scholar]

[128] Shalgar, S., & Tamborra, I. 2024, Phys. Rev. D, 109, 103011 [NASA ADS] [CrossRef] [Google Scholar]

[129] Shibata, M., Kiuchi, K., Sekiguchi, Y.-i., & Suwa, Y. 2011, Prog. Theor. Phys., 125, 1255 [CrossRef] [Google Scholar]

[130] Smith, M. C., Sijacki, D., & Shen, S. 2018, MNRAS, 478, 302 [NASA ADS] [CrossRef] [Google Scholar]

[131] Sotani, H., & Takiwaki, T. 2016, Phys. Rev. D, 94, 044043 [NASA ADS] [CrossRef] [Google Scholar]

[132] Sotani, H., Kuroda, T., Takiwaki, T., & Kotake, K. 2017, Phys. Rev. D, 96, 063005 [NASA ADS] [CrossRef] [Google Scholar]

[133] Sotani, H., Kuroda, T., Takiwaki, T., & Kotake, K. 2019, Phys. Rev. D, 99, 123024 [CrossRef] [Google Scholar]

[134] Sotani, H., Takiwaki, T., & Togashi, H. 2021, Phys. Rev. D, 104, 123009 [NASA ADS] [CrossRef] [Google Scholar]

[135] Steiner, A. W., Hempel, M., & Fischer, T. 2013, ApJ, 774, 17 [NASA ADS] [CrossRef] [Google Scholar]

[136] Sukhbold, T., Woosley, S. E., & Heger, A. 2018, ApJ, 860, 93 [NASA ADS] [CrossRef] [Google Scholar]

[137] Suleiman, L., Oertel, M., & Mancini, M. 2023, Phys. Rev. C, 108, 035803 [NASA ADS] [CrossRef] [Google Scholar]

[138] Sullivan, C., O’Connor, E., Zegers, R. G. T., Grubb, T., & Austin, S. M. 2016, ApJ, 816, 44 [NASA ADS] [CrossRef] [Google Scholar]

[139] Sumiyoshi, K., Yamada, S., Suzuki, H., et al. 2005, ApJ, 629, 922 [NASA ADS] [CrossRef] [Google Scholar]

[140] Summa, A., Janka, H. T., Melson, T., & Marek, A. 2018, ApJ, 852, 28 [NASA ADS] [CrossRef] [Google Scholar]

[141] Takiwaki, T., & Kotake, K. 2018, MNRAS, 475, L91 [NASA ADS] [CrossRef] [Google Scholar]

[142] Takiwaki, T., Kotake, K., & Suwa, Y. 2014, ApJ, 786, 83 [NASA ADS] [CrossRef] [Google Scholar]

[143] Thielemann, F.-K., Isern, J., Perego, A., & von Ballmoos, P. 2018, Space Sci. Rev., 214, 62 [CrossRef] [Google Scholar]

[144] Thorne, K. S. 1980, Rev. Mod. Phys., 52, 299 [NASA ADS] [CrossRef] [Google Scholar]

[145] Torres-Forné, A., Cerdá-Durán, P., Passamonti, A., & Font, J. A. 2018, MNRAS, 474, 5272 [CrossRef] [Google Scholar]

[146] Torres-Forné, A., Cerdá-Durán, P., Obergaulinger, M., Müller, B., & Font, J. A. 2019a, Phys. Rev. Lett., 123, 051102; Erratum: 2021, 127, 239901 [CrossRef] [Google Scholar]

[147] Torres-Forné, A., Cerdá-Durán, P., Passamonti, A., Obergaulinger, M., & Font, J. A. 2019b, MNRAS, 482, 3967 [CrossRef] [Google Scholar]

[148] Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9 [NASA ADS] [CrossRef] [Google Scholar]

[149] Vartanyan, D., Burrows, A., Radice, D., Skinner, M. A., & Dolence, J. 2018, MNRAS, 477, 3091 [NASA ADS] [CrossRef] [Google Scholar]

[150] Vartanyan, D., Burrows, A., Radice, D., Skinner, A. M., & Dolence, J. 2019, MNRAS, 482, 351 [NASA ADS] [CrossRef] [Google Scholar]

[151] Vartanyan, D., Coleman, M. S. B., & Burrows, A. 2022, MNRAS, 510, 4689 [NASA ADS] [CrossRef] [Google Scholar]

[152] Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]

[153] Wang, T., & Burrows, A. 2024, ApJ, 962, 71 [NASA ADS] [CrossRef] [Google Scholar]

[154] Westernacher-Schneider, J. R., O’Connor, E., O’Sullivan, E., et al. 2019, Phys. Rev. D, 100, 123009 [NASA ADS] [CrossRef] [Google Scholar]

[155] Wolfe, N. E., Fröhlich, C., Miller, J. M., Torres-Forné, A., & Cerdá-Durán, P. 2023, ApJ, 954, 161 [NASA ADS] [CrossRef] [Google Scholar]

[156] Woosley, S. E., & Heger, A. 2007, Phys. Rep., 442, 269 [Google Scholar]

[157] Xiong, Z., Wu, M.-R., Abbar, S., et al. 2023a, Phys. Rev. D, 108, 063003 [NASA ADS] [CrossRef] [Google Scholar]

[158] Xiong, Z., Wu, M.-R., Martínez-Pinedo, G., et al. 2023b, Phys. Rev. D, 107, 083016 [CrossRef] [Google Scholar]

[159] Yasin, H., Schäfer, S., Arcones, A., & Schwenk, A. 2020, Phys. Rev. Lett., 124, 092701 [NASA ADS] [CrossRef] [Google Scholar]

[160] Zappa, F., Bernuzzi, S., Radice, D., & Perego, A. 2023, MNRAS, 520, 1481 [Google Scholar]