11institutetext: The Oskar Klein Centre, Department of Astronomy
Stockholm University, AlbaNova, SE-106 91 Stockholm, Sweden 22institutetext:
Department of Physics and Astronomy, Michigan State University,
East Lansing, MI 48824, USA 33institutetext: Department of Computational
Mathematics, Science, and Engineering, Michigan State University,
East Lansing, MI 48824, USA 44institutetext: Facility for Rare Isotope Beams,
Michigan State University,
East Lansing, MI 48824, USA

# Grey Two-moment Neutrino Transport: Comprehensive Tests and Improvements for
Supernova Simulations

Haakon Andresen 11<EMAIL_ADDRESS>Evan P. O’Connor 11 Oliver
Eggenberger Andersen 11 and Sean M. Couch 223344

###### Abstract

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

Methods. We implemented the framework of Foucart et al. (2016) 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 grey 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 grey 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 per cent
discrepancies in the average energy and luminosity of heavy-lepton neutrinos.
Simulations using the grey transport are around four times faster than those
using energy-dependent transport.

###### Key Words.:

supernovae: general – Neutrinos – Radiative transfer – Hydrodynamics –
Gravitational waves

##  1 Introduction

Core-collapse supernovae are important for several facets of astrophysics,
they play 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 so far
by LIGO-Virgo-KAGRA collaboration (Abbott et al., 2016b, a, c, 2017b, 2017a,
2017c, 2017d, 2019, 2020c, 2020b, 2020a, 2020d, 2021a, 2021b, 2021b, 2024,
2023). 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 modelling 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; Burrows et al., 2019; Vartanyan et al.,
2019; Burrows et al., 2020; Vartanyan et al., 2022). 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; Martin et al., 2021;
Johns, 2023; Capozzi et al., 2022; Xiong et al., 2023b; Richers & Sen, 2022;
Dedin Neto et al., 2023; Liu et al., 2023; Xiong et al., 2023a; Cornelius et
al., 2023; Shalgar & Tamborra, 2023; Ehring et al., 2023; Dedin Neto et al.,
2023; Cornelius et al., 2023; Shalgar & Tamborra, 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 aforementioned aspects. The input
physics underlying core-collapse supernovae are complex and subject to active
research. Consequently, the underpinning of supernova simulations comes with
inherent uncertainty.

On 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 their 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, yet 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; Summa et al., 2018; Vartanyan et al., 2018; Kuroda et al., 2018;
Burrows et al., 2020; Wang & Burrows, 2023; Kuroda et al., 2022), flux limited
diffusion (Bruenn et al., 2020; Rahman et al., 2019, 2022), fast multi-group
transport (Müller & Janka, 2015), and the isotropic diffusion source
approximation (Takiwaki et al., 2014; 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 (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 their numerical cost and
complexity, neither full Boltzmann methods nor Monte Carlo methods have been
widely adopted in simulations (see Mezzacappa et al. (2020) for a detailed
review). Energy-integrated, or grey, schemes are often employed in simulations
of binary neutron star mergers (Foucart et al., 2015, 2016; Fujibayashi et
al., 2020; Radice et al., 2022; Hayashi et al., 2022; Zappa et al., 2023;
Kiuchi et al., 2023a; Fujibayashi et al., 2023; Curtis et al., 2023; Radice &
Bernuzzi, 2023; Kiuchi et al., 2023b; Schianchi et al., 2023; Foucart, 2023).
While grey transport methods have been 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).

Grey 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 grey 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 grey scheme to a more sophisticated method, we improve
several aspects of the original grey transport and achieve good overall
agreement between the two approaches. We start by implementing the energy-
integrated neutrino transport suggested by Foucart et al. (2015) and Foucart
et al. (2016) into FLASH. We then test its accuracy and suggest specific
changes based on our test results.

This paper is organised as follows: We first briefly describe the current
version of FLASH in section 2. In section 3 we summarise 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 section 4. In section 5 we present the changes
we propose to the original transport scheme. Towards the end of section 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 6. In section 7 we present the
gravitational waves (GWs) and neutrino signals from a simulation utilising our
grey scheme and compare the results to signals from a corresponding simulation
with energy-dependent transport. In section 8 we discuss the performance of
the grey transport, compared to the energy-dependent transport. Finally, in
section 9 we summarise 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 utilise a modified version of the FLASH
framework (version 4) (Fryxell et al., 2000), which has been 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 of 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 are generated with the neutrino opacity library NuLib (O’Connor,
2015; Sullivan et al., 2016). Heavy-lepton neutrinos are mainly produced in
the PNS through pair-production processes, the rates of which 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. An effective absorption is calculated through Kirchhoff’s law
from the effective emissivity. We refer the reader 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),

$\displaystyle p^{\alpha}\bigg{[}\frac{\partial f_{\nu}}{\partial
x^{\alpha}}-\Gamma^{\beta}_{\alpha\gamma}p^{\gamma}\frac{\partial
f_{\nu}}{\partial p^{\beta}}\bigg{]}$ $\displaystyle=\bigg{[}\frac{\partial
f_{\nu}}{\partial\tau}\bigg{]}_{\mathrm{coll}}.$ (1)

Here $f=f(x^{\alpha},p^{\alpha})$ is the distribution function of the
neutrinos, with $x^{\alpha}=(t,x^{i})$ and $p^{\alpha}$ representing the four-
position and the four-momentum of the neutrinos, respectively.
$\Gamma^{\beta}_{\alpha\gamma}$ 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 $\nu$ represents the different neutrino species, in this work we
evolve three neutrino species: electron neutrinos ($\nu_{e}$), electron anti-
neutrinos ($\bar{\nu}_{e}$), and a third species ($\nu_{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 grey, moments.

###  3.1 Evolution equations

We summarise the method proposed by Foucart et al. (2016), but we refer the
reader 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
section 5, we will 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,

$\displaystyle\partial_{t}E+\partial_{i}\big{[}\alpha F^{i}\big{]}+\alpha
F^{i}\partial_{i}\Phi=\alpha\bigg{(}W[\eta-\kappa_{a}J]-\kappa_{H}H^{t}\bigg{)},$
(2) $\displaystyle\partial_{t}F^{i}+\partial_{j}[\alpha P^{ij}]+\alpha
E\partial_{i}\Phi=\alpha\bigg{(}W\big{[}\eta-\kappa_{a}J\big{]}v^{i}-\kappa_{H}H^{i}\bigg{)},$
(3) $\displaystyle\partial_{t}N+\partial_{j}\big{[}\alpha
F_{N}^{j}\big{]}=\alpha\bigg{(}{\eta}_{N}-\frac{{\kappa}_{N}J{N}}{W({E}-{F}_{i}v^{i})}\bigg{)}.$
(4)

Here $\eta$ represents the energy-integrated emissivity, $\eta_{N}$ is the
energy-integrated number emissivity, and $\kappa_{N}$ is the energy-averaged
number absorption. $\kappa_{H}$ denotes the sum of the energy-averaged
absorption ($\kappa_{a}$) and energy-averaged scattering opacities
($\kappa_{s}$), $\kappa_{H}=\kappa_{a}+\kappa_{s}$. The energy-averaged
opacities are detailed in section 3.3. The fluid velocity is represented by
$v^{i}$, $W$ is the Lorentz factor, the gravitational potential is denoted
with $\Phi$, and the lapse is $\alpha=\exp{(\Phi)}$. The terms $J$ and
$H^{\mu}$ denote the energy density and momentum density in the fluid rest
frame, respectively, and are related to the laboratory variables as follows:

$\displaystyle J$ $\displaystyle=W^{2}[E-2F^{i}v_{i}+v^{i}v^{j}P_{ij}],$ (5)
$\displaystyle H^{t}$
$\displaystyle=W^{3}[-(E-F^{i}v_{i})v^{2}+F^{i}v_{i}-v_{i}v_{j}P^{ij}],$ (6)
$\displaystyle H^{i}$
$\displaystyle=W^{3}[-(E-F^{j}v_{j})v^{i}+F^{i}-v_{j}P^{ij}].$ (7)

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

In Eq. (4), $F^{j}_{N}$ is the number momentum density that we assume takes
the form of

$F^{j}_{N}=\frac{JWv^{j}}{\langle\varepsilon\rangle}+\frac{H^{j}}{\langle\varepsilon^{F}\rangle},$
(8)

where $\langle\varepsilon\rangle$ is the average neutrino energy and
${\langle\varepsilon^{F}\rangle}$ is the average energy of the neutrino flux
(see section 3.4).

###  3.2 Neutrino spectrum

The grey scheme of Foucart et al. (2016) is underpinned by the assumption that
the neutrinos follow a Fermi-Dirac distribution

$\displaystyle
f_{\nu}(\varepsilon)=\frac{1}{1+\exp((\varepsilon-\mu_{\nu})/T_{\nu})},$ (9)

where $\varepsilon$ represents the neutrino energy, $\mu_{\nu}$ is the
chemical potential of the neutrinos and is the $T_{\nu}$ the neutrino
temperature. With the corresponding black body intensity

$\displaystyle
B_{\nu}(\varepsilon)\propto\frac{\varepsilon^{3}}{1+\exp((\varepsilon-\mu_{\nu})/T_{\nu})}.$
(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

$\langle\varepsilon\rangle=W\frac{E-F_{i}v^{i}}{N}.$ (11)

For a black body 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

$\langle\varepsilon\rangle=\frac{F_{3}(\eta_{\nu})}{F_{2}{(\eta_{\nu}})}T_{\nu},$
(12)

where $\eta_{\nu}=\mu_{\nu}/T$ (not to be confused with the emissivity or
number emissivity, $\eta$ or $\eta_{N}$). The Fermi integrals are defined as
follows 111These integrals can be numerically challenging and we are required
to evaluate a large number of them for each time step. We follow the accurate
and fast approach of Fukushima (2015) to compute the integrals.:

$\displaystyle
F_{k}(\eta_{\nu})=\int_{0}^{\infty}\frac{x^{k}}{1+e^{x-\eta_{\nu}}}\mathrm{d}x.$
(13)

We can then find the neutrino temperature by combining equations Eq. (11) and
Eq. (12). In our implementation of the grey scheme, we follow Foucart et al.
(2015) and chose the following expression for the chemical potential

$\displaystyle\mu_{\nu}=\mu_{\nu}^{eq}(1-e^{-\tau}),$ (14)

where $\mu_{\nu}^{eq}$ is the chemical potential of the neutrinos assuming
equilibrium with the fluid and $\tau$ is the optical depth of the neutrinos.

###  3.3 Neutrino Opacities

The energy-integrated neutrino-matter source terms in Eqs. 2, 3, and 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

$\displaystyle\eta=\int\eta(\rho,Y_{e},T,\varepsilon)\mathrm{d}\varepsilon.$
(15)

In the above equation, $\rho$, $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,

$\displaystyle\kappa_{a}=\frac{\int\kappa_{a}(\varepsilon)J(\varepsilon)\mathrm{d}\varepsilon}{\int
J(\varepsilon)\mathrm{d}\varepsilon},$ (16)

and

$\displaystyle\kappa_{H}=\frac{\int\kappa_{H}(\varepsilon)H^{r}(\varepsilon)\mathrm{d}\varepsilon}{\int
H^{r}(\varepsilon)\mathrm{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 black body 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\propto B_{\nu}$, and define

$\displaystyle\kappa_{a/H}^{eq}=\frac{\int
B_{\nu}(T,\mu_{\nu},\varepsilon)\kappa_{a/H}(\rho,Y_{e},T,\varepsilon)\mathrm{d}\varepsilon}{\int
B_{\nu}(T,\mu_{\nu},\varepsilon)\mathrm{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

$\displaystyle\kappa_{a/H}=\kappa_{a/H}^{eq}\frac{T_{\nu}^{2}}{T^{2}}.$ (19)

The correction factor is based on the fact that the opacities typically go as
energy squared, in other words, $\kappa\sim\varepsilon^{2}$. Defining the
opacities in terms of Eq. (19) ensures that

$\kappa_{a}=\eta J$ (20)

in the optically thick regions where the fluid is in equilibrium with the
fluid ($T_{\nu}=T$) and consequently guarantees detailed balance.

For the number emissivity and number absorption Foucart et al. (2016) chose

$\kappa_{N}=\kappa_{a}\frac{\eta_{N}}{\eta}\frac{F_{3}(\eta_{\nu})T}{F_{2}(\eta_{\nu})},$
(21)

with $\kappa_{a}$ from Eq. (19) and

$\displaystyle\eta_{N}=\int\eta(\rho,Y_{e},T,\varepsilon)\varepsilon^{-1}\mathrm{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 grey counterparts are calculated using NuLib (O’Connor,
2015; Sullivan et al., 2016). A Python script to generate energy-integrated
NuLib tables is available at https://github.com/evanoconnor/NuLib

###  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 grey 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 $\varepsilon^{2}$, one
expects (Rosswog & Liebendoerfer, 2003)

$\displaystyle\langle\varepsilon^{F}\rangle=\frac{F_{1}(\eta_{\nu})}{F_{0}(\eta_{\nu})}T_{\nu}.$
(23)

Based on the expected flux-weighted average energy, Foucart et al. (2016)
introduced two new scalars and wrote down a simple equation for modelling the
ratio of the flux-weighted average energy and the average neutrino energy:

$\displaystyle\frac{\langle\varepsilon^{F}\rangle}{\langle\varepsilon\rangle}=\frac{F_{3}(\eta_{\nu})F_{0}(\eta_{\nu})-s^{F}(F_{3}(\eta_{\nu})F_{0}(\eta_{\nu})-F_{2}(\eta_{\nu})F_{1}(\eta_{\nu}))}{F_{3}(\eta_{\nu})F_{0}(\eta_{\nu})-s^{C}(F_{3}(\eta_{\nu})F_{0}(\eta_{\nu})-F_{2}(\eta_{\nu})F_{1}(\eta_{\nu}))}.$
(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 towards zero in the optically thick region, where the neutrinos
are in equilibrium with the fluid, and towards $s^{F}$ in the optically thin
regions, where the neutrinos are streaming, one possible choice is

$\displaystyle
s^{C}=\frac{Ns^{C}_{0}+\alpha\mathcal{F}s^{F}dt}{N+\alpha\mathcal{F}dt+\bar{\eta}_{N}\alpha
dt}.$ (25)

In Eq. (25), $s_{0}^{C}$ is the value at the beginning of the time step and
$\mathcal{F}$ is an estimate for the neutrino number flux. The scalar $s_{F}$
can be defined as follows

$\displaystyle s^{F}=\frac{s^{C}+\tau}{1+\tau},$ (26)

The optical depth can be estimated based on the flux factor
($\xi=\sqrt{H^{\mu}H_{\mu}}/J$) and a free parameter $\beta$ as follows:

$\xi=\frac{1}{1+\beta\tau}.$ (27)

The approximate neutrino number flux was set to

$\mathcal{F}=\xi
N\bigg{(}\frac{F_{3}(\eta_{\nu})F_{0}(\eta_{\nu})-s^{F}(F_{3}(\eta_{\nu})F_{0}(\eta_{\nu})-F_{2}(\eta_{\nu})F_{1}(\eta_{\nu}))}{F_{2}(\eta_{\nu})F_{1}(\eta_{\nu})}\bigg{)}^{2}.$
(28)

##  4 Comparison of Foucart et al. (2016) to Energy-dependent Transport

To evaluate the accuracy of the original grey scheme, we perform a set of
tests in one dimension (1D). The tests were all performed with $\beta=6$ (in
Eq. (27)). We start 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.

Table 1: Neutrino luminosities and average neutrino energy in hydrostatic simulations of the post bounce phase of a 30 solar mass progenitor. Each row represents a different neutrino transport scheme. The top row shows the results in the simulation with energy-dependent transport, and the following rows represent different iterations of the grey neutrino-transport scheme. Luminosities are denoted by $L$ and average energies by $\langle\varepsilon\rangle$ with a subscript to denote the three different neutrino species. | $L_{\nu_{e}}\,[10^{51}\mathrm{erg/s}]$ | $L_{\bar{\nu}_{e}}\,[10^{51}\mathrm{erg/s}]$ | $L_{\nu_{x}}\,[10^{51}\mathrm{erg/s}]$ | $\langle\varepsilon_{\nu_{e}}\rangle\,[\mathrm{MeV}]$ | $\langle\varepsilon_{\bar{\nu}_{e}}\rangle\,[\mathrm{MeV}]$ | $\langle\varepsilon_{\nu_{x}}\rangle\,[\mathrm{MeV}]$
---|---|---|---|---|---|---
ED | 91.62 | 28.84 | 129.29 | 10.34 | 12.11 | 19.02
Foucart et al. (2016) | 82.55 | 23.00 | 82.48 | 8.72 | 10.90 | 22.83
Grey | 91.51 | 29.56 | 117.19 | 9.11 | 11.12 | 25.93

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 grey transport (which we will return to below). The original grey
transport shows relatively good agreement with the energy-dependent neutrino
transport. The average energies of all neutrino species agree to within 20 per
cent in the two methods. The luminosities of the electron-type neutrinos agree
within $\sim 10$ per cent, the electron anti-neutrino luminosities agree to
$\sim 20$ per cent, and the heavy-lepton neutrino luminosities differ by
almost 40 per cent. The luminosities in the simulation using the original grey
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 500 ms 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,2, and 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 dull-coloured lines.
The solid purple line shows the results from the simulation with our
implementation of the grey transport from Foucart et al. (2016).

Figure 1: Neutrino luminosities for the grey simulation compared to the
simulations presented in O’Connor et al. (2018) (dull coloured lines). The
purple lines indicate the results from the simulation based on the original
grey 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
(dotted-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 cut-out in the left panel shows a zoom-in
of the peak centred around $\sim 0.1$ s after bounce.

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 towards the end of the simulation.

Figure 2: Average neutrino energies for the grey simulation (bright lines)
compared to the simulations presented in O’Connor et al. (2018) (dull coloured
lines). The purple lines indicate the results from the simulation based on the
original grey 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 (dotted-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.

The transport of Foucart et al. (2016) leads to average neutrino energies for
the electron neutrinos and electron anti-neutrinos that are within $\sim 20$
per cent 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 around 60 per cent larger compared to the typical value seen in
the other simulations (see the right panel of Fig. 2).

Figure 3: Shock radius (solid lines) and PNS radius (dashed lines) for the
grey simulation (purple lines) compared to the simulations presented in
O’Connor et al. (2018) (dull colored lines). The purple lines indicate the
results from the simulation based on the original grey 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^{11}$ g/cm3.

From Fig. 3, we see that both the PNS radius and the shock radius are
significantly larger in the grey 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 grey neutrino transport.

##  5 Improved grey 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 $\kappa_{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 non-linear, it can be challenging to construct comparisons that
isolate a single component in dynamical simulations. We, therefore, base 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 grey 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 perform the following
evaluation: Using the spectral properties of the radiation fields in the
simulation with energy-dependent neutrino transport, we compute

$\displaystyle\kappa_{a}^{ed}$ $\displaystyle=\frac{\int
J(\varepsilon)\kappa_{a}(\rho,Y_{e},T,\varepsilon)\mathrm{d}\varepsilon}{\int
J(\varepsilon)\mathrm{d}\varepsilon},$ (29)

and

$\displaystyle\kappa_{H}^{ed}$
$\displaystyle=\frac{\int\big{[}\kappa_{a}(\rho,Y_{e},T,\varepsilon)+\kappa_{s}(\rho,Y_{e},T,\varepsilon)\big{]}H^{r}(\varepsilon)\mathrm{d}\varepsilon}{\int
H^{r}(\varepsilon)\mathrm{d}\varepsilon},$ (30)

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

Figure 4: Energy-averaged absorption opacity as a function of radius for a
hydrostatic evolution of a 30 solar mass progenitor using the SFHo EOS. Dashed
red lines show the opacities calculated with the neutrino spectrum from the
energy-dependent neutrino transport (Eq. (29)). The solid blue lines show the
opacities calculated according to Eq. (19). The dashed green 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 Figs. 4 and 5, we show the opacity values for $\kappa_{a}$ and
$\kappa_{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 develop in this work, which we will return to later in this
section. The figures show $\kappa_{a}$ and $\kappa_{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 section 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 grey 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 black-body
assumption consistently yields higher values of $\kappa_{H}$ than the red line
based on the full neutrino spectrum.

Figure 5: The 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. 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.

The correction factor in Eq. 19 is an extrapolation away from equilibrium
between the fluid and neutrinos. The anchor point $\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 employ an interpolation method with pre-computed values of the
opacities defined as follows

$\displaystyle\kappa_{a/H}(T_{\nu})=\frac{\int
B_{\nu}(T_{\nu},\mu_{\nu},\varepsilon)\kappa_{a/H}(\rho,Y_{e},T,\varepsilon)\mathrm{d}\varepsilon}{\int
B_{\nu}(T_{\nu},\mu_{\nu},\varepsilon)\mathrm{d}\varepsilon}.$ (31)

We chose a set of 12 linearly spaced neutrino temperatures from 0.1 MeV to 16
MeV. 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_{\nu}>16\,$MeV was
beneficial from a numerical standpoint.

The interpolation method shows 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 found that the interpolation approach solved an
issue we had noticed in our implementation of the original grey 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 neutrinosphere 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$ will lead 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 222The 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.. Since we do not
have information regarding the energy dependence of the momentum flux in our
energy-integrated scheme, we multiply $\kappa_{H}$ with $0.6$ in the core to
mimic this effect.

Implementing the correction to $\kappa_{H}$ led to improved agreement between
the PNS and neutrino properties of the grey scheme and the full energy-
dependent neutrino transport. However, the grey 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 $\kappa_{a}$ and $\kappa_{H}$ as follows:

$\displaystyle\bar{\kappa}_{a}$
$\displaystyle=\bigg{(}\cos^{2}\big{(}{\frac{\pi\xi}{2}}\big{)}+0.95\sin^{2}\big{(}{\frac{\pi\xi}{2}}\big{)}\bigg{)}\kappa_{a}(T_{\nu}),$
(32)

and

$\displaystyle\bar{\kappa}_{H}$
$\displaystyle=\bigg{(}0.6\cos^{2}\big{(}{\frac{\pi\xi}{2}}\big{)}+0.95\sin^{2}\big{(}{\frac{\pi\xi}{2}}\big{)}\bigg{)}\kappa_{H}(T_{\nu}).$
(33)

Here $\kappa_{a}(T_{\nu})$ and $\kappa_{H}(T_{\nu})$ refers to the
interpolated values of the $\kappa_{a}$ and $\kappa_{H}$, respectively. We
note that the factor of $0.6$ does not appear in absorption opacity, which
shows up in the term $W[\eta-\kappa_{a}J]$ 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{\kappa}_{a}$ 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 Fig.
4 and Fig. 5, respectively. For the heavy-lepton neutrinos, we see good
agreement between $\bar{\kappa}_{a}$ and ${\kappa}_{a}^{ed}$. The absorption
opacities for the electron neutrinos and electron anti-neutrinos agree
relatively well between all three methods, but $\bar{\kappa}_{a}$ is slightly
closer to ${\kappa}_{a}^{ed}$ than ${\kappa}_{a}^{eq}$. For $\kappa_{H}$, we
see that $\bar{\kappa}_{H}$ agrees better with the value from the energy-
dependent calculation than ${\kappa}_{H}^{eq}$. However, there is still a some
discrepancy between $\bar{\kappa}_{H}$ and ${\kappa}_{H}^{ed}$, for all
neutrino species.

Our tests showed that getting better results for $\kappa_{H}$ is possible by
choosing a different function and different coefficients in Eq. (33). However,
the results shown in Fig. 4 and Fig. 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, Eq. (32) and Eq. (33) were chosen to produce the best
overall agreement in all our tests (see 5.3, and 6).

###  5.2 Updated Flux-weighted Average Neutrino Energy

To improve the machinery described in section 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

$\displaystyle\frac{F_{1}(\eta_{\nu})}{F_{0}(\eta_{\nu})}T_{\nu}<\langle\varepsilon_{F}\rangle<\frac{F_{3}(\eta_{\nu})}{F_{2}(\eta_{\nu})}T_{\nu}.$
(34)

Results from the energy-dependent transport indicate that the flux-weighted
average energy matches well with

$\displaystyle\langle\varepsilon_{F}\rangle=\frac{F_{2}(\eta_{\nu})}{F_{1}(\eta_{\nu})}T_{\nu},$
(35)

especially for the heavy-lepton neutrinos for which the second term in Eq. (8)
is particularly important. We, therefore, change the order of the Fermi
integrals in Eq. (24) and Eq. (28) to better match the results from the
energy-dependent transport. The updated expressions are

$\mathcal{F}=\xi
N\bigg{(}\frac{F_{3}(\eta_{\nu})F_{1}(\eta_{\nu})-s^{F}(F_{3}(\eta_{\nu})F_{1}(\eta_{\nu})-F_{2}(\eta_{\nu})F_{2}(\eta_{\nu}))}{F_{2}(\eta_{\nu})F_{2}(\eta_{\nu})}\bigg{)}^{2},$
(36)

and

$\displaystyle\frac{\langle\varepsilon^{F}\rangle}{\langle\varepsilon\rangle}=\frac{F_{3}(\eta_{\nu})F_{1}(\eta_{\nu})-s^{F}(F_{3}(\eta_{\nu})F_{1}(\eta_{\nu})-F_{2}(\eta_{\nu})F_{2}(\eta_{\nu}))}{F_{3}(\eta_{\nu})F_{1}(\eta_{\nu})-s^{C}(F_{3}(\eta_{\nu})F_{1}(\eta_{\nu})-F_{2}(\eta_{\nu})F_{2}(\eta_{\nu}))}.$
(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
$\langle\varepsilon^{F}\rangle=\langle\varepsilon\rangle$ (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

$\displaystyle
s^{C}=\frac{Ns^{C}_{0}+\alpha\mathcal{F}s^{F}dt}{N(1+\tau)+\alpha\mathcal{F}dt+\bar{\eta}_{N}\alpha
dt}.$ (38)

The only difference between Eq. (38) and Eq. (25) is the term $N(1+\tau)$ 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 $\tau$ as follows

$\displaystyle\tau(r)=\int_{r}^{\infty}\langle\bar{\kappa}_{H}\rangle_{\Omega}^{-1}\mathrm{d}r,$
(39)

where $\langle\bar{\kappa}_{H}\rangle_{\Omega}$ represents the angular average
of $\bar{\kappa}_{H}$ at any given radius. We use $\bar{\kappa}_{H}$ to
calculate the optical depth since we are attempting to estimate the average
energy of the neutrino flux. Directly calculating $\tau$ is advantageous since
it completely does away with the $\beta$ parameter of Foucart et al. (2016), a
parameter that influences the original scheme’s results. Foucart et al. (2016)
reported a $5$-$10$ per cent variation in the neutrino luminosites and average
energies for reasonable values of $\beta$ ($\in(0,8]$). We found a similar
dependence on the $\beta$ parameter, as Foucart et al. (2016), in our
implementation of their scheme. While calculating the optical depth directly
improves the result of our grey scheme, relying on the approximate value from
Eq. (27) is not detrimental to the accuracy of the transport.

The accuracy of our grey scheme depends on the two free parameters we
introduced in Eq. (32) and Eq. (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 five per cent.

###  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 per cent. The average energy of the heavy-
lepton neutrinos agrees within 30 per cent. 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 grey scheme.

The black curves in Figs. 1,2, and 3 show the results using the updated grey
scheme. The new grey 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; Andresen
et al., 2019; Morozova et al., 2018; Westernacher-Schneider et al., 2019;
Sotani et al., 2021; Eggenberger Andersen et al., 2021; Vartanyan et al.,
2022; Mori et al., 2023).

##  6 Two-dimensional tests

We extended our testing of the grey scheme by conducting 2D axisymmetric core-
collapse supernova simulations, using both the fully energy-dependent neutrino
transport and the grey transport.

Figure 6: The 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. The dashed red line
corresponds to simulations with fully energy-dependent transport, and the
green lines indicates simulations using the grey scheme of Foucart et al.
(2016), and blue lines show results from simulations using the grey 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.

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:

  1. 1.

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 labelled 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.

  2. 2.

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
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 the reader 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 grey 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 grey
simulations were mapped from a spherically symmetric profile 20 ms after
bounce. It is necessary to start the grey 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 grey transport and simulations using energy-dependent neutrino
transport for a wide range of progenitor and for several different EOSs (see
Fig. 6 and Fig. 7.)

Figure 7: The 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, green lines indicate simulations using the grey
scheme of Foucart et al. (2016), and blue lines show results from simulations
using the grey 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.

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^{11}$g/cm3) 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 grey 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 utilising energy-
dependent transport.

From the right-hand side panels 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 grey 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
$\sim 33\,$km, $\sim 35.5\,$km, and $\sim 38\,$km in the simulations using
energy-dependent transport, our grey implementation, and the grey transport
according to Foucart et al. (2016), respectively. For the s14 progenitor,
again at 400 ms post bounce, we find a PNS radius of $\sim 34$ km using
energy-dependent transport, $\sim 36.5$ km using our grey transport and $\sim
39$ km using the original grey 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 grey 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 grey transport. The
grey scheme is more conducive to shock revival, and shock expansion tends to
set in sooner in simulations using the grey transport than in simulations
using energy-dependent transport. The situation is notably reversed for the
original grey scheme; the original grey 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\times 10^{10}\,$ g/cm3 and the
entropy is greater than 6 $k_{\mathrm{B}}$ per baryon ($k_{\mathrm{B}}$
denotes Boltzmann constant). In Fig. 8, blue curves indicate results from
simulations with our grey 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 grey scheme injects too much energy into
the fluid at early times, the heating rates are to high, with the notable
exception of model s12.

Figure 8: The heating rate as a function of time for models s11, s12, s13,
s14, and s15. Blue curves indicate results from simulations with our grey
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. Each row correspond to
one model, the model names are shown in the bottom left corner of each panel.

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 grey transport and the energy-dependent
transport. Especially since the evolution of the grey 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 behaviour 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 grey scheme indicate that our
scheme could be improved in the future.

Table 2: Overview of the model names for our second set of 2D simulations. The table shows the relationship between the model name the effective mass parameter, and the isoscalar incompressibility modulus parameter. The effective mass is given in units of the free neutron mass $m_{n}$. Model name | $m^{\star}$ | $K_{\mathrm{sat}}$
---|---|---
| [$m_{n}$] | [MeV baryon-1]
s20_0.55 | 0.55 | 230
s20_0.75_200 | 0.75 | 200
s20_0.75 | 0.75 | 230
s20_0.75_260 | 0.75 | 260
s20_0.95 | 0.95 | 230

We compare the simulations from the second set in Fig. 7. We plot the average
shock radius and the PNS radius for each realisation 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 grey 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
utilising 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 grey scheme
consistently results in larger PNS radii than the energy-dependent transport.
Furthermore, we see that the updated grey transport again leads to larger PNS
radii than the energy-dependent transport, from approximately $\sim 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. (2016) leads
to later shock revival and our grey scheme leads to earlier explosions
compared to the energy-dependent transport. However, shock revival sets in
within 50 ms in the simulation using our grey 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
grey transport. It is absent in the simulation that incorporates the modified
grey scheme. Similar behaviour 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 grey transport produces qualitatively different behaviour for
all three models. On the other hand, the transport of Foucart et al. (2016)
yields better agreement across the three EOS variations. The 2D tests show
that the grey transport shifts the outcome towards explosion for progenitors
in the boundary region between explodability and non-explodability.
Consequently, our findings warrant caution in using grey transport to predict
the outcome of individual progenitors precisely. However, our grey 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.

##  7 Gravitational Waves and Neutrino Emission

In this section, we present the neutrino luminosity, average energy, and GW
signal for the energy-dependent and grey simulations based on the s12
progenitor. We do not show results for the grey implementation of Foucart et
al. (2016). 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_{\times}$, but only one component is non-zero 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

$\displaystyle h_{+}=\frac{3}{2}\frac{G}{c^{4}D}\ddot{Q}_{33}$ (40)

Here is the $c$ the speed of light and $G$ is Newton’s constant.
$\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{\mathrm{d}}{\mathrm{d}t}\bigg{[}\int\,\rho
zv_{z}\mathrm{d}^{3}x\bigg{]}.$ (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
$\rho$ 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 $(-\infty,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.

Figure 9: Gravitational waves 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 grey 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 have been normalized by a common factor.

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 grey 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 $\sim$ 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 grey 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 towards lower frequencies in the
simulation employing grey 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 grey 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; Torres-Forné et al., 2018; Sotani et al.,
2019; Morozova et al., 2018; Radice et al., 2019; Torres-Forné et al., 2019b,
a; Sotani et al., 2021; 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 towards
lower frequencies in the grey transport simulation.

Figure 10: Neutrino luminosities (top panel) and average energies (bottom)
from two simulations of the s12 progenitor. Red lines indicate results from
the simulation with energy-dependent neutrino transport. Blue lines show
results from the simulation performed with our grey neutrino transport. Solid
lines indicate electron neutrinos, dashed lines show electron anti-neutrinos,
and dotted-dashed lines correspond to heavy-lepton neutrinos.

Fig. 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 grey 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 grey neutrino
transport induces notable differences for the heavy-lepton neutrinos compared
to the energy-dependent transport. Initially, the grey transport produces
approximately 10 MeV higher energies, with around 30 per cent lower
luminosities than those observed in the energy-dependent transport simulation.
At late times, the luminosity agreement improves to about 15-20 per cent.

The $\sim 10\,$MeV offset in average-energy of the neutrinos heavy-lepton,
between the simulations with grey and energy-integrated transport, is present
for all our models. Furthermore, all of our grey 2D models show the same 15-30
per cent 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 grey 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.

The differences we observe in the neutrino emission are consistent with the
results in 1D test simulations, see section 5.

##  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 grey 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 hydrodynamic calculation. For
comparison, in the energy-dependent scheme, a radiation step is $\sim 15$
times more expensive than a hydrodynamic step.

##  9 Summary and Conclusion

In this work, we have 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 per cent 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 grey
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 grey scheme. We implemented the following changes to the approach
proposed by Foucart et al. (2016):

  1. 1.

Updated the definition of energy-averaged neutrino opacities. Concretely:

     * •

Foucart et al. (2016) computes 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_{\nu}/T)^{2}$. We assume that the neutrinos obey a black body
spectrum (which was also assumed by Foucart et al. (2016)) and compute the
opacities for a set of predefined neutrino temperatures and create an
interpolation grid in $T_{\nu}$.

     * •

After interpolating the opacities, we apply 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 $(\kappa_{a}+\kappa_{s})H^{i}$ with 0.6 in optically thick
regions.

  2. 2.

Adjusted the procedure for estimating the average energy of the neutrino flux.
Specifically, we made the following changes

     * •

Changed the order of the Fermi integrals in the equations modelling the ratio
of the average energy of the neutrino flux and the average neutrino energy.
For consistency, we also changed the updated the order of the Fermi integrals
that appear in the estimate for the neutrino number flux.

     * •

Modified the scalar $s^{C}$ by introducing an $N\tau$ term in the denominator.

     * •

Directly computed the optical depth, eliminating reliance on the $\beta$
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 grey 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
while the neutrino temperature was high. Interpolating the opacities, instead
of extrapolating, fixed this problem. Not only does our approach reproduce the
results of the energy-dependent transport fairly well, simulations utilising
the grey 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
$\big{(}\kappa_{s}+\kappa_{a}\big{)}H^{i}$, 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 tests could be solved by reducing the opacities by 5 per cent 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 travelling
through an optically thick screen will be reduced from
$\frac{F_{3}(\eta_{\nu})}{F_{2}(\eta_{\nu})}T_{\nu}$ to
$\frac{F_{1}(\eta_{\nu})}{F_{0}(\eta_{\nu})}T_{\nu}$. Empirically, we found
that the energy reduction is closer to
$\frac{F_{2}(\eta_{\nu})}{F_{1}(\eta_{\nu})}T_{\nu}$. 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, which 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 grey scheme. Overall, we found
that for a range of progenitors and for different EOSs the grey 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 grey transport and
simulations with multi-group M1 transport. However, the grey scheme is
slightly more conducive to shock revival and we conclude that it is likely to
push simulations residing in the boundary region between successful and failed
supernovae over the edge towards explosion.

We performed two sets of 2D test simulations. In the first set we tested the
accuracy of the grey 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 grey 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 $\sim 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 grey 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 grey scheme will favour explosions in cases close to the boundary between
failed and successful supernova explosions. Again, we found that the original
grey transport hampers shock revival and leads to extended PNSs.

One motivation for implementing a grey 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 grey transport produces GW and neutrino
signals that agree well with those of the energy-dependent transport. In
section 7, we compared the predicted GW and neutrino signals for one
particular model using the grey transport with the corresponding model with
energy-dependent transport. The GW signals agreed well, but we observed a
small shift towards lower frequencies in the grey model, which we attributed
to the larger PNS radius observed in the grey model (a tendency that is
systematic to the grey 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 grey scheme, was expected.

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