# Neutrino flavor mixing with moments

McKenzie Myers  Department of Physics, North Carolina State University,
Raleigh, NC 27695, USA  Theo Cooper  Department of Physics, North Carolina
State University, Raleigh, NC 27695, USA  MacKenzie Warren  AI Research Group,
5x5 Technologies, St Petersburg, FL 33701, USA  Jim Kneller
<EMAIL_ADDRESS>Department of Physics, North Carolina State University,
Raleigh, NC 27695, USA  Gail McLaughlin  Department of Physics, North Carolina
State University, Raleigh, NC 27695, USA  Sherwood Richers  NSF Astronomy &
Astrophysics Postdoctoral Fellow  Department of Physics, University of
California Berkeley, Berkeley, CA 94720  Evan Grohs  Department of Physics,
North Carolina State University, Raleigh, NC 27695, USA  Carla Fröhlich
Department of Physics, North Carolina State University, Raleigh, NC 27695, USA

###### Abstract

The successful transition from core-collapse supernova simulations using
classical neutrino transport to simulations using quantum neutrino transport
will require the development of methods for calculating neutrino flavor
transformations that mitigate the computational expense. One potential
approach is the use of angular moments of the neutrino field, which has the
added appeal that there already exist simulation codes which make use of
moments for classical neutrino transport. Evolution equations for quantum
moments based on the quantum kinetic equations can be straightforwardly
generalized from the evolution of classical moments based on the Boltzmann
equation. We present an efficient implementation of neutrino transformation
using quantum angular moments in the free streaming, spherically symmetric
bulb model. We compare the results against analytic solutions and the results
from more exact multi-angle neutrino flavor evolution calculations. We find
that our moment-based methods employing scalar closures predict, with good
accuracy, the onset of collective flavor transformations seen in the multi-
angle results. However in some situations they overestimate the coherence of
neutrinos traveling along different trajectories. More sophisticated quantum
closures may improve the agreement between the inexpensive moment-based
methods and the multi-angle approach.

##  I Introduction

The importance of neutrinos to the core-collapse supernova (CCSN) paradigm has
been recognized since the earliest simulations of the explosions by Colgate
and White [1]. Even though many details have changed over the years since that
pioneering study, neutrinos are still thought to be the driver of the
explosions of stars with initial masses $\gtrsim 10\,M_{\odot}$ that have
reached the end of their nuclear burning lifetimes (see [2, 3, 4, 5, 6, 7] for
reviews).

Although a complete theoretical accounting of the explosion dynamics remains
elusive, state of the art numerical simulations performed in three dimensions
appear to be converging toward exploding solutions [8, 9, 10, 11, 12, 13, 14,
15, 16]. However, there are still many uncertainties underlying these models.
The structure of the progenitor stars is presently poorly understood and is
complicated by variations in stellar mass, metallicity, rotation history,
magnetic fields, and binary interactions. Another source of uncertainty in the
details of the explosion arises from the unknown equation of state describing
matter above nuclear densities. Variations in the equation of state have been
shown to lead to significant differences in the outcomes of simulations (e.g.,
[17, 18, 19]). Furthermore, there is yet work to be done to ensure that
multidimensional simulations are performed with sufficient resolution and that
the approximations employed in various methods (especially regarding neutrino
transport) do not significantly affect the solution [20, 21], though there
currently seems to be more agreement than disagreement between full supernova
simulation codes [22, 23].

While the recent progress in supernova simulations is impressive, the
simulations neglect the quantum nature of the neutrino and its ability to
change flavor (though see [24, 25] for some attempts at effective treatments).
Until recently, this omission was not thought to be important for the dynamics
of the explosion (see [26, 27] for recent reviews). Flavor transformation
calculations based on post-processing results from simulations that employ
classical transport reveal a number of neutrino flavor transformation
phenomena. As in the Sun, neutrinos undergo complete flavor transformation at
a Mikheyev-Smirnov-Wolfenstein (MSW) resonance [28, 29, 30] due to the
combination of differing neutrino masses and a potential from neutrinos
interacting with the background matter. However, this flavor transformation
occurs well outside of any supernova engine or nucleosynthesis region and is
mostly important for understanding future neutrino detections at Earth, e.g.
[31, 32, 33]. The neutrino self-interaction potential due to neutrinos
interacting with other nearby neutrinos makes the flavor transformation
nonlinear. The so-called collective neutrino oscillations, which result from
the combined effects of the self-interaction potential and the differing
neutrino masses, are also thought to occur too far out to affect the supernova
engine in a significant way [34, 35, 24], unless beyond the standard model
interactions are included, e.g. [36].

But more recently, the so-called fast flavor instability was shown to have the
potential to drive significant and rapid flavor transformations well inside
the shock if the neutrino distribution fulfils the criterion of an electron
lepton number crossing [37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 47].
Conditions fulfilling this instability criterion have been shown to exist both
inside [49, 50, 51, 52] and outside [53] the shock front in simulations that
use classical transport. Recent consideration of effects of small-scale
turbulence could make the fast flavor instability even more ubiquitous than
the electron lepton number crossing criterion predicts [54, 55, 56]. Thus
these more recent studies now strongly indicate that flavor transformation may
occur in regions of the supernova where there could be feedback into the
hydrodynamics. While simulations that adopt effective treatments of flavor
transformation within the framework of classical neutrino transport may be
able to capture such physics in an approximate way, direct calculations of
quantum neutrino transport are needed to understand the effects of microscopic
instabilities from first principles.

There are several challenges to be overcome in implementing and simulating
quantum neutrino transport. One unavoidable challenge is a huge increase in
the dynamical range of lengthscales that need to be resolved. At the present
time, the spatial resolution of the state-of-the-art simulations approaches
lengthscales of order $\sim~{}100\;{\rm m}$ for a simulation covering a domain
size of order $\sim 1000\;{\rm km}$. However, the lengthscale of neutrino
flavor oscillations in matter with a density of $\sim 10^{12}\;{\rm g/cm^{3}}$
is of order $\sim 1\;{\rm\mu m}$, which means one would need a simulation with
a dynamical range of at least twelve orders of magnitude. Multi-angle
collective flavor transformation calculations also require hundreds to
thousands of angle bins in order for the results to converge, far more than
the $\sim 10$ angle bins often used in simulations of classical transport [57,
58, 21, 59]. The results from multi-angle calculations are often seen to
exhibit substantial changes seen as the number of angle bins used is changed.
The number of energy bins used in multi-angle calculations is also usually of
order a few hundred, an order of magnitude larger than the $\sim 20$ energy
bins typically used in classical simulations. However, in contrast to the
number of angle bins, the convergence of the multi-angle results with the the
number of energy bins is often observed to be much smoother [58] and the large
number of energy bins used is driven by the desire for sufficient resolution
of the spectrum. Thus the computational expense of even a 1D, spherically
symmetric supernova simulation using quantum transport is expected to be many
orders of magnitude greater than for a simulation using classical transport.
One must seriously consider alternative approaches that mitigate this expense
if quantum supernova simulations are to be feasible.

One approach, which we consider here, is to use angular moments of the quantum
neutrino distribution i.e. quantum moments. There are good motivations for
considering this approach. First, the coupling of the neutrinos to the rest of
the fluid is most simply expressed via the moments so computing the moments
directly is more efficient than computing the flavor evolution of neutrinos
traveling along different trajectories and then integrating. Second, neutrino
transport based on angular moments of the classical neutrino radiation field
are already used in many state of the art supernova simulation codes due to
their computational efficiency (e.g ., [60, 61, 62, 63, 64, 65]). Modifying
such codes to include the equations describing the evolution of the quantum
moments will rely on techniques developed both in the neutrino oscillation
literature and in the neutrino transport literature. Vlasenko  _et al._ [66]
and Volpe [67] developed the equations from first principles, while Blaschke
and Cirigliano [68] expressed the collision integral in detail. Investigations
of the feasibility of the quantum moment approach have already been taken.
Strack and Burrows [69] (see also Duan and Shalgar [70]) outlined a moment-
based formalism for the QKEs and Richers  _et al._ [71] fleshed out the form
of the collision integral using interaction rates in the form commonly used in
the core-collapse supernova simulation literature. Richers  _et al._ [71] also
developed a code to simulate the QKEs under the assumption of isotropy and
homogeneity, and Johns  _et al._ [72] used a moment method to analyze the
presence of fast flavor instability in parameterized cases. When comparing
moment-based approaches to more-sophisticated transport, disagreement emerges
from the need to truncate the tower of moment evolution equations at some
level yielding one or two fewer equations than the number of unknown moments.
To solve the equations one must propose an algorithm to estimate the unknown
moment(s) given the evolved moments. If this algorithm were perfect, the
evolved moments are guaranteed to exactly match those extracted from a full
Boltzmann calculation. Truncating the tower of moment equations at different
levels does not necessarily lead to more accurate results. The accuracy of the
results is determined by the accuracy of the algorithm for finding the un-
evolved moment(s).

The goal of this paper is to present results from a new code that solves the
quantum moment evolution equations for a supernova neutrino bulb model and
allows us to explore these issues. In Section II we present the quantum
kinetic moment equations we use, introduce the one-moment and two-moment
schemes and explain the closures we adopt for them. In Section III we
investigate how well moment methods are able to produce MSW and collective
oscillations by comparing their results to multi-angle calculations. Finally,
we conclude in Section IV and discuss the issues our study uncovers.

##  II Flavor mixing with moments

Classical radiation transport is described in general by the Boltzmann
equation, but the six-dimensional distribution function is exceedingly
computationally expensive to evolve. A moment approach to radiation transport
involves taking angular moments of the full Boltzmann equation and evolving a
subset of these moments rather than directly discretizing the neutrino
distribution in angle. The result is a system of coupled differential
equations for a handful of four-dimensional (three spatial and one energy
dimension) variables, rather than a differential equation for a six-
dimensional distribution function. The angular moments which are evolved or
prescribed in a two-moment method in flat spacetime are the differential
energy density $E$, (energy) flux vector $\vec{F}$, and pressure tensor
$\tensor{P}$, defined as

$\displaystyle E(t,\vec{r},q)$
$\displaystyle=\frac{1}{4\pi}\,\left(\frac{q}{2\pi\hbar c}\right)^{3}\int
d\Omega_{p}\,f(t,\vec{r},\vec{p})$ (1) $\displaystyle\vec{F}(t,\vec{r},q)$
$\displaystyle=\frac{1}{4\pi}\,\left(\frac{q}{2\pi\hbar c}\right)^{3}\int
d\Omega_{p}\,\hat{p}\,f(t,\vec{r},\vec{p})$ (2)
$\displaystyle\tensor{P}(t,\vec{r},q)$
$\displaystyle=\frac{1}{4\pi}\,\left(\frac{q}{2\pi\hbar c}\right)^{3}\int
d\Omega_{p}\,\hat{p}\otimes\hat{p}\,f(t,\vec{r},\vec{p})\,\,,$ (3)

where $q=|\vec{p}|\,c$. Note that we omit the speed of light from the
definition of the flux moment in order that the moments all have the same
units. The factor of $1/4\pi$ in each definition of the moments is a units
choice made to divide out the $4\pi$ from the solid angle integral when
tracking the evolution of these variables. This means, the total neutrino flux
$\vec{\cal{F}}$ is $4\pi$ times the integral of the spectral flux

${\vec{\cal{F}}}=4\pi\,\int_{0}^{\infty}{\vec{F}}\,dq\,\,.$ (4)

The moments for the antineutrinos are defined similarly and shall be denoted
with overbars i.e. $\bar{E}$, $\vec{\bar{F}}$ and $\tensor{\bar{P}}$. In
classical transport $f$ and $\bar{f}$ are distribution functions, taking on
values $f\in[0,1]$ and describing the occupation number per unit of phase
space. To generalize so as to permit flavor mixing, we refine $f$ and
$\bar{f}$ to become matrices in flavor space. Thus, the moments also become
matrices and the energy density, for example, is now (assuming two neutrino
flavors) [69]

$E(t,\vec{r},q)\rightarrow\begin{pmatrix}E^{(ee)}(t,\vec{r},q)&E^{(ex)}(t,\vec{r},q)\\\
E^{(xe)}(t,\vec{r},q)&E^{(xx)}(t,\vec{r},q)\end{pmatrix}\,\,.$ (5)

The other moments have a similar structure. We shall refer to moments that are
matrices in flavor space as ‘quantum moments’ in order to distinguish them
from the scalar ‘classical’ definition of moments, and indicate flavor
elements of the quantum moments with superscripts in parentheses to avoid
confusions with exponents. $E^{(ee)}$ and $E^{(xx)}$ are real, positive
definite quantities that represent the differential energy densities in the
electron and representative heavy lepton (“$x$”) flavor neutrinos. The off-
diagonal components $E^{(ex)}$ and $E^{(xe)}$ represent the flavor overlap and
are complex quantities. The luminosity matrices of neutrinos and antineutrinos
are calculated from the radial component $F_{r}$ and $\bar{F}_{r}$ of the flux
vectors for the neutrinos and antineutrinos respectively, and are

$\displaystyle L$ $\displaystyle=\left(4\pi r\right)^{2}c\int F_{r}\,dq,$ (6)
$\displaystyle\bar{L}$ $\displaystyle=\left(4\pi
r\right)^{2}c\int\bar{F}_{r}\,dq\,\,.$ (7)

Similarly to the quantum moments for the energy density, flux and pressure,
the diagonal elements of the luminosity matrices are the luminosities of each
neutrino flavor.

The evolution equations for the moments are found by following the procedure
outlined in Strack and Burrows [69] and Zhang and Burrows [73] (although the
collision terms therein need to be modified as in [68, 71]). There are an
infinite number of angular moments and corresponding evolution equations (e.g.
[74]). Under the assumption of flat spacetime and spherical symmetry, the
evolution equations for the lowest-order moments (energy density and energy
flux) take the form

$\displaystyle\frac{1}{c}\frac{\partial E}{\partial t}+\frac{\partial
F_{r}}{\partial r}+\frac{2F_{r}}{r}$
$\displaystyle=-\frac{i}{\hbar\,c}\left[H_{V}+H_{M}+H_{E},E\right]-\frac{i}{\hbar\,c}\left[H_{F},F_{r}\right]+C_{E}$
(8) $\displaystyle\frac{1}{c}\frac{\partial\bar{E}}{\partial
t}+\frac{\partial\bar{F}_{r}}{\partial r}+\frac{2\bar{F}_{r}}{r}$
$\displaystyle=-\frac{i}{\hbar\,c}\left[H_{V}-H_{M}-H_{E}^{*},\bar{E}\right]+\frac{i}{\hbar\,c}\left[H_{F}^{*},\bar{F}_{r}\right]+\bar{C}_{E}$
(9) $\displaystyle\frac{1}{c}\,\frac{\partial F_{r}}{\partial
t}+\frac{\partial P_{rr}}{\partial r}+\frac{3\,P_{rr}-E}{r}$
$\displaystyle=-\frac{i}{\hbar\,c}\left[H_{V}+H_{M}+H_{E},F_{r}\right]-\frac{i}{\hbar\,c}\left[H_{F},P_{rr}\right]+C_{F}$
(10) $\displaystyle\frac{1}{c}\,\frac{\partial\bar{F}_{r}}{\partial
t}+\frac{\partial\bar{P}_{rr}}{\partial r}+\frac{3\,\bar{P}_{rr}-\bar{E}}{r}$
$\displaystyle=-\frac{i}{\hbar\,c}\left[H_{V}-H_{M}-H_{E}^{*},\bar{F}_{r}\right]+\frac{i}{\hbar\,c}\left[H_{F}^{*},\bar{P}_{rr}\right]+\bar{C}_{F}$
(11)

where we have defined

$H_{E}=\phantom{-}4\pi\,\sqrt{2}\,G_{F}\,\int\frac{dq}{q}\left(E(q)-\bar{E}^{*}(q)\right)$
(12)
$H_{F}=-4\pi\,\sqrt{2}\,G_{F}\,\int\frac{dq}{q}\left(F_{r}(q)-\bar{F}_{r}^{*}(q)\right)$
(13)

and ∗ indicates the complex conjugate. Together $H_{E}$ and $H_{F}$ are the
self-interaction Hamiltonian $H_{SI}=H_{E}+H_{F}$. The four terms $C_{E}$,
$\bar{C}_{E}$, $C_{F}$ and $\bar{C}_{F}$ are the ‘collision’ terms and, in the
context of flavor mixing, are also matrices [68, 71]. We neglect the collision
terms in the rest of this work. For two flavor mixing, the vacuum Hamiltonian
has the form

$H_{V}=\frac{\Delta
m^{2}_{12}\,c^{4}}{4\,q}\,\left[\,\sin(2\,\theta_{12})\,\sigma_{1}-\cos(2\,\theta_{12})\,\sigma_{3}\,\right]\,\,,$
(14)

where $\Delta m_{12}^{2}=m_{2}^{2}-m_{1}^{2}$ is the splitting between the
neutrino masses and $\theta_{12}$ is the neutrino vacuum mixing angle. The
matter Hamiltonian has the form

$H_{M}=\frac{\sqrt{2}G_{F}n_{e}}{2}\sigma_{3}~{},$ (15)

where $n_{e}$ is the number density of electrons. In both cases, $\sigma_{i}$
is the $i$’th Pauli matrix and we have removed terms that contribute only to
the trace of the Hamiltonian, as they do not affect the flavor evolution.

In order to solve for the evolution of the moments we need to truncate the
tower of equations. In classical moment transport this is typically done at
the first or second level i.e. Eqns. 8 and 9, or 8 to 11. But this truncation
introduces a hurdle: inspection of the truncated tower of equations reveals
that they contain, in general, one more moment than the number of equations
allows us to solve for. Thus to solve the truncated tower we need to close the
evolution equations by divining a relationship between the moments that are
evolved and the moments which are not. A perspicacious choice for the closure
will result in the agreement of moments with higher resolution methods.
Conversely, choosing an approximate closure will only yield approximate
solutions. Furthermore, truncating the tower of equations at a higher level
does not necessarily lead to more accurate solutions for the moments that are
solved, unless the structure of the higher moments is known better than that
of the lower moments. Rather, evolving more moments simply provides more
information with which to construct closures, though increased accuracy of the
solutions is not guaranteed just because we use a closure with more inputs.

###  II.1 Neutrino Bulb Model

In order to test the ability of moment-based schemes to reproduce the flavor
transformation from more sophisticated transport calculations, we need to
compare results from the two approaches. Unfortunately there are not many test
problems which allow such a comparison. The test cases we consider in this
work focus upon the inhomogeneous environment around a spherical source of
neutrinos i.e. a neutrino bulb. Since this environment has been studied
extensively due to its similarity to core-collapse supernovae, we can take
advantage of results from other codes which have been developed for this
scenario. In what follows we compare our results with the steady-state ‘multi-
angle’ BULB model calculations using the SQA code, which is briefly described
in Appendix A.

For the moment-based calculations we shall consider truncations at the first
and second level of the equation tower i.e. a scheme where we solve for $E$
and $\bar{E}$, or $E$, $\bar{E}$, $F_{r}$ and ${\bar{F}}_{r}$. For each we
need a closure. In previous studies of classical moment transport, closures
have been supplied from analytic physical approximations or ad-hoc
prescriptions (see [75, 76] for summaries) or characteristic methods (e.g.,
[77, 78]). There is a wealth of physics buried in the choice of closure and
any results will depend on this choice. We leave exploration of the
sensitivity to the closure to future work and use here two geometrically
motivated examples.

Figure 1: The bulb model: at a distance $r$ from the center of a opaque
neutrino bulb of radius $R_{\nu}$, the neutrino has traveled a distance
$\lambda$ and its trajectory makes an angle $\theta$ with the radial
direction.

The geometry of the scenario we are considering is shown in Figure 1. At a
given radial location $r$ above the source, the neutrinos are confined to
propagate within a cone around the radial direction. Neutrinos which are
propagating along a trajectory which makes an angle $\theta$ with the radial
direction have traveled a distance

$\lambda=r\cos\theta-\sqrt{R_{\nu}^{2}-r^{2}\,\sin^{2}\theta}$ (16)

from the neutrinosphere $R_{\nu}$. The half opening-angle of the cone
containing all the neutrino trajectories passing through the point of interest
is $\theta_{\mathrm{max}}$ and from the geometry this angle is found to be

$\cos\theta_{\mathrm{max}}=\sqrt{1-(R_{\nu}/r)^{2}}.$ (17)

We will use this to self-consistently inform our choice of closures.

####  II.1.1 One-moment closure

The first closure is a relationship between the energy and the radial
component of the flux for use in one moment calculations. In classical flux-
limited diffusion, the system is similarly closed so as to directly evolve
only one moment, though the specific closure is different from that employed
here. When the neutrinos are emitted half-isotropically from a spherical bulb
and do not oscillate, the integrals defining the scalar energy density and
radial flux (or more generally, the traces thereof) are related analytically.
In this case we find the flux is related to the energy density via

$F_{r}=\frac{(1-\cos^{2}\theta_{\mathrm{max}})}{2\,(1-\cos\theta_{\mathrm{max}})}E~{}.$
(18)

Using Equation 17, this can be used to define a closure as

$F_{r}=\frac{E}{2}\left(1+\sqrt{1-\left(\frac{R_{\nu}}{r}\right)^{2}}\right).$
(19)

In what follows we will refer to this relationship as the “one-moment”
closure. Substituting the one-moment closure of Eq. 19 into Eq. 8 gives an
equation that is now only a function of the two energy densities. Using this
closure, the form of the self-interaction Hamiltonian matches that in the
single-angle approximation of Dasgupta  _et al._ [79], although we stress that
the advection terms in the one-moment approximation differ from those in the
single-angle approximation.

####  II.1.2 Two moment closure

In the two moment approximation we simultaneously solve for the evolution of
the first two moments as is done in M1 neutrino transport methods common in
CCSN modeling. These moments, the energy density $E$ and flux $F_{j}$, must be
related to the next highest moment, the pressure tensor $P_{ij}$. In the
neutrino bulb model, in the absence of oscillations, we find that the
$rr$-component of the pressure tensor can be related to the energy density via

$P_{rr}=\frac{(1-\cos^{3}\theta_{\mathrm{max}})}{3\,(1-\cos\theta_{\mathrm{max}})}E~{}.$
(20)

Using Equation 17, this becomes

$P_{rr}=\left(2-\left(\frac{R_{\nu}}{r}\right)^{2}+\sqrt{1-\left(\frac{R_{\nu}}{r}\right)^{2}}\right)\frac{E}{3}.$
(21)

In what follows we shall refer to this equation as the “two-moment“ closure.

####  II.1.3 Inner Boundary Condition

To construct the initial values for the moments at the neutrinosphere we adopt
a steady state distribution $f$ which is pure diagonal in the flavor basis.
The diagonal entries of the distribution matrix are of the form
$f^{(aa)}(R_{\nu},q,\theta)\propto Q_{a}(q)\,\Theta_{a}(\theta)$ with $Q_{a}$
a function describing the energy spectrum for flavor $a$ and $\Theta_{a}$ a
function describing the angular spectrum for flavor $a$. For simplicity, the
energy spectra are taken to be Fermi Dirac distributions

$Q_{a}(q)=\frac{1}{1+\exp{\left(q/(k_{B}\,T_{a})-\eta_{a}\right)}}$ (22)

where $T_{a}$ is the temperature and $\eta_{a}$ is the chemical potential
divided by $k_{B}T_{a}$. For both the moment-based and multi-angle
calculations we adopt a uniform energy grid from $1\,\mathrm{MeV}$ to
$60\,\mathrm{MeV}$ with 591 energy bins corresponding to an energy resolution
of $100\;{\rm keV}$. Repeating all calculations using a resolution of
$50\;{\rm keV}$ indicates the numerical error is less than 0.1%. We adopt the
parametric form for the initial angular distributions introduced by Mirizzi
and Serpico [80]

$\Theta_{a}(\theta)=\cos^{\beta_{a}}\theta.$ (23)

The case of half isotropic emission corresponds to $\beta_{a}=0$. Using these
distributions we write

$f^{(aa)}(R_{\nu},q,\theta)=A_{a}\,(\beta_{a}+2)\,Q_{a}(q)\,\Theta_{a}(\theta)$
(24)

and imposing the requirement that the luminosity matrix $L$ be given by
equation (6) at the neutrinosphere $R_{\nu}$, we find the constant $A_{a}$ to
be

$A_{a}=\frac{\pi}{G_{2}(\eta)}\left(\frac{c}{R_{\nu}}\right)^{2}\,\left(\frac{\hbar}{k_{B}T_{a}}\right)^{3}\,\frac{L^{(aa)}}{\langle
q\rangle_{a}}\,\,.$ (25)

The function $G_{2}(\eta)$ is the complete Fermi-Dirac integral defined as

$G_{j}(x)=\int_{0}^{\infty}\frac{t^{j}}{1+\exp{\left(t-x\right)}}dt$

and the quantity $\langle q\rangle_{a}$ is the mean energy given by

$\langle q\rangle_{a}=\frac{\int F^{(aa)}(R_{\nu},q)\,dq}{\int
F^{(aa)}(R_{\nu},q)/q\,dq}~{}.$ (26)

With the distribution matrix now defined, the energy density and flux moments
at radius r (assuming no oscillations) are computed to be

$\displaystyle E(r,q)$
$\displaystyle=\frac{R_{\nu}^{2}}{2\,r^{2}}\,\left(\frac{q}{2\pi\hbar
c}\right)^{3}\begin{pmatrix}A_{e}\,Q_{e}(q)\,_{2}F_{1}\left(1,\frac{1}{2},2+\frac{\beta_{e}}{2},\frac{R_{\nu}^{2}}{r^{2}}\right)&0\\\
0&A_{x}\,Q_{x}(q)\,_{2}F_{1}\left(1,\frac{1}{2},2+\frac{\beta_{x}}{2},\frac{R_{\nu}^{2}}{r^{2}}\right)\end{pmatrix}$
(27) $\displaystyle F_{r}(r,q)$
$\displaystyle=\frac{R_{\nu}^{2}}{2\,r^{2}}\,\left(\frac{q}{2\pi\hbar
c}\right)^{3}\begin{pmatrix}A_{e}\,Q_{e}(q)&0\\\
0&A_{x}\,Q_{x}(q)\end{pmatrix}$ (28) $\displaystyle P_{rr}(r,q)$
$\displaystyle=\frac{R_{\nu}^{2}}{2\,r^{2}}\,\left(\frac{q}{2\pi\hbar
c}\right)^{3}\begin{pmatrix}A_{e}\,Q_{e}(q)\,_{2}F_{1}\left(1,-\frac{1}{2},2+\frac{\beta_{e}}{2},\frac{R_{\nu}^{2}}{r^{2}}\right)&0\\\
0&A_{x}\,Q_{x}(q)\,_{2}F_{1}\left(1,-\frac{1}{2},2+\frac{\beta_{x}}{2},\frac{R_{\nu}^{2}}{r^{2}}\right)\end{pmatrix}$
(29)

where ${}_{2}F_{1}(a,b,c,z)$ is the ordinary hypergeometric function. Similar
expressions give the initial conditions for the antineutrino moments.

###  II.2 Numerical Method

To solve the quantum moment evolution equations we have developed a new code
that computes the steady-state solution based on an inner boundary condition,
given a specified matter density and electron fraction throughout the
computational domain. We omit the collision terms in the moment evolution
equations in order to focus our attention upon the oscillation physics. The
code solves the equations using an explicit midpoint (second order Runge-
Kutta) integrator. The radial step is adaptive based on the three frequencies
associated with the various flavor mixing terms in the Hamiltonian:

$\displaystyle\omega_{\mathrm{V}}$ $\displaystyle=\frac{\Delta
m_{12}^{2}\,c^{4}}{2\,\hbar\,q}$ (30) $\displaystyle\omega_{\mathrm{M}}$
$\displaystyle=\sqrt{2}G_{F}\,n_{e}/\hbar$ $\displaystyle\omega_{\mathrm{SI}}$
$\displaystyle=\sqrt{||H_{SI}||^{2}_{F}-\frac{1}{2}\,\left[Tr(H_{SI})\right]^{2}}\,/\hbar$

and $||H_{SI}||^{2}_{F}$ is the Frobenius norm of the self-interaction
Hamiltonian, and $Tr(H_{SI})$ its trace. In the absence of the collision
terms, the quantities
$||F_{r}||^{2}_{F}-\frac{1}{2}\,\left[Tr(F_{r})\right]^{2}$ and
$||\bar{F}_{r}||^{2}_{F}-\frac{1}{2}\,\left[Tr(\bar{F}_{r})\right]^{2}$ are
conserved, i.e. independent of the radius $r$, for every energy. Our code
enforces the conservation of these quantities to an error tolerance set to
$0.1\%$ per step. If the fractional change in the size of these quantities
exceeds this bound, the time increment is halved and the step repeated. In
what follows we present the results from several test problems of the code.

##  III Testing Moment-Based Methods

If moment-based approaches to flavor oscillations are to be a feasible
alternative to calculations based on discrete ordinates or other more exact
approaches, the results must agree well with analytic predictions if they
exist and/or the results from less approximate numerical approaches.

###  III.1 Constant Electron Density

Our first test case is flavor mixing of neutrinos emitted from a hard sphere
of radius $R_{\nu}$ in a background of constant electron density. The flavor
evolution for a single neutrino is well known for this scenario. The
probability $P_{T}(r,\theta)$ that a neutrino initially in a particular flavor
transitions to a neutrino of the opposite flavor after traveling a distance
$\lambda$ is given by a sinusoid with fixed amplitude and wavelength.
Specifically,

$P_{T}(r,\theta)=\sin^{2}\left(2\theta_{\mathrm{MSW}}\right)\sin^{2}\left(\frac{\omega_{\mathrm{MSW}}\,\lambda}{\hbar\,c}\right)~{}.$
(31)

In this equation the effective matter mixing angle $\theta_{\mathrm{MSW}}$ and
effective frequency $\omega_{\mathrm{MSW}}$, are

$\displaystyle\sin^{2}(2\,\theta_{\mathrm{MSW}})$
$\displaystyle=\frac{\sin^{2}(2\,\theta_{12})}{\sin^{2}(2\,\theta_{12})+C^{2}}$
(32) $\displaystyle\omega_{\mathrm{MSW}}$
$\displaystyle=\frac{\omega_{\mathrm{V}}}{2}\sqrt{\sin^{2}(2\,\theta_{12})+C^{2}}~{},$

where

$C=\cos(2\theta_{12})\mp\frac{\omega_{\mathrm{M}}}{\omega_{\mathrm{V}}}~{}.$
(33)

The negative sign is used for neutrinos and the positive sign for
antineutrinos.

Using equation (31), one can compute the angular moments of an ensemble of
neutrinos by integrating the single neutrino solution over the solid angle
with the appropriate $\cos\theta$ weight and taking into account the different
path lengths $\lambda$ from the emission point to the radial point of
interest. However one does not need to do this calculation to predict what one
should observe in the solutions. Due to the differing path lengths, the
neutrinos will not all have the same phase of the flavor oscillations (i.e.,
they will be not all be coherent). The incomplete coherence will inevitably
smear the oscillations of the moments compared to the flavor oscillations of a
single neutrino. The amount of coherence between the neutrinos on different
trajectories is a function of the distance from the source. When the spread in
path lengths to a given point is much smaller than the wavelength of the
oscillations, all the neutrinos are very close to being coherent, but as one
moves away from the source the spread in path lengths grows and the amount of
coherence drops. For a neutrinosphere of radius $R_{\nu}$, the path length
difference between neutrinos emitted radially and those emitted at a tangent
to the neutrinosphere asymptotes to $R_{\nu}$. Thus the amount of coherence
should also asymptote to a level determined by the ratio of the oscillation
wavelength to $R_{\nu}$.

In order to test our code’s ability to reproduce these predictions we solve
the two moment evolution equations and separately perform a multi-angle
calculation using SQA. The spatial grid extends to $40\,\mathrm{km}$, we use
artificial vacuum mixing parameters of $\Delta m^{2}_{12}=6.9\times
10^{-4}\,\mathrm{eV}^{2}$ and $\theta_{12}=0.28818$ and consider monoenergetic
neutrinos with an energy of $q=1\,\mathrm{MeV}$. The elements of the energy
density moments are set arbitrarily according to the hierarchy
$E^{(ee)}(R_{\nu})=10\,\bar{E}^{(ee)}(R_{\nu})=100\,E^{(xx)}(R_{\nu})=100\,\bar{E}^{(xx)}(R_{\nu})$.
The matter density is set to $\rho=8\times 10^{3}\,\mathrm{g\,cm}^{-3}$ and
the electron fraction is $Y_{e}=0.5$ in order to put the $1\,\mathrm{MeV}$
neutrinos close to the MSW resonance.

Figure 2: Flavor oscillation test including vacuum and matter contributions to
the Hamiltonian. Here we use $\Delta m^{2}_{12}=6.9\times
10^{-4}\,\mathrm{eV}^{2}$ and $\theta_{12}=0.28818$ in a background density
profile of $\rho=8000\,\mathrm{g/cm}^{3}$ and $Y_{e}=0.5$. The top panel shows
the transition probability of electron neutrinos to $x$-flavor neutrinos while
the bottom is the same for the antineutrinos. Solid lines are the results from
the moment code and dashed lines are the results from the SQA using 90,001
angular bins.

Figure 2 shows the oscillation probabilities as a function of the radius $r$
for two angular distribution parameters of ${\beta_{a}}=0$ (dashed red) and
${\beta_{a}}=3$ (dashed purple). The electron neutrino and antineutrino
transition probabilities are defined from the flux moment to be

$\displaystyle P_{\nu_{e}\rightarrow\nu_{x}}(r)$ $\displaystyle=$
$\displaystyle\frac{r^{2}\,F_{r}^{(xx)}(r)-R_{\nu}^{2}\,F_{r}^{(xx)}(R_{\nu})}{R_{\nu}^{2}\,[F_{r}^{(ee)}(R_{\nu})-F_{r}^{(xx)}(R_{\nu})]}$
(34) $\displaystyle{\bar{P}}_{\bar{\nu}_{e}\rightarrow\bar{\nu}_{x}}(r)$
$\displaystyle=$
$\displaystyle\frac{r^{2}\,\bar{F}_{r}^{(ee)}(r)-R_{\nu}^{2}\,\bar{F}_{r}^{(xx)}(R_{\nu})}{R_{\nu}^{2}\,[\bar{F}_{r}^{(ee)}(R_{\nu})-\bar{F}_{r}^{(xx)}(R_{\nu})]}.$
(35)

The top panel shows the flavor evolution as a function of radius for the
neutrinos and the bottom panel for the antineutrinos. Although the neutrinos
are on resonance and undergo nearly complete flavor oscillations,
antineutrinos are off-resonance and only undergo minor flavor transformation.
A comparison of the SQA results for the half-isotropic case and a separate,
analytic integration of the survival probability, i.e.,

$P_{\nu_{e}\rightarrow\nu_{x}}(r)=\frac{2\,r^{2}}{R_{\nu}^{2}}\,\int^{1}_{\cos\theta_{max}}P_{T}(r,\theta)\,\cos\theta\,d(\cos\theta)\,\,,$
(36)

are visually indistinguishable.

The SQA results (dashed) display the decoherence effect of the neutrinos that
have traveled along different trajectories leading to the reduction in the
amplitude of the oscillations with increasing radius. Close to the
neutrinosphere, the amplitudes of the flavor oscillations predicted from the
moment code (solid) are similar in magnitude to the ${\beta_{a}}=0$ SQA
results (dashed purple). However, as the neutrinos move away, the moment code
maintains larger flavor oscillations than exhibited in the multi-angle
results.

Although the moment code overestimates the coherence of the neutrinos
traveling along different trajectories, it does capture some phase effects. As
the distribution becomes more forward-peaked, the average phase advances more
slowly due to smaller average path length to a given radial point. This can be
seen by comparing the SQA results for ${\beta_{a}}=3$ (forward-peaked) and
${\beta_{a}}=0$ (semi-isotropic). The moment results, which have initial
conditions corresponding to ${\beta_{a}}=0$, show a more slowly evolving phase
than the ${\beta_{a}}=3$ results, just as the SQA ${\beta_{a}}=0$ results do.
However, the agreement does not last for long and a significant phase
difference between the the moment and ${\beta_{a}}=0$ oscillations exists at
larger radius. Of course, in the limit of a perfectly forward-peaked
distribution (not shown), the moment equations become increasingly accurate
and agree increasingly well with multi-angle data.

Figure 3: The ratio of the ‘ee’ element of luminosity matrix relative to its
initial value, versus radius for the test in Section III.2.1 and Table 1. The
top panel is for neutrinos and the bottom panel is for antineutrinos. The one-
moment results are difficult to discern due to being largely obscured by the
two-moment results.

###  III.2 The Steady State Self Interaction Problem In Spherical Symmetry

Our next three tests are for the case of collective neutrino oscillations in
an inhomogeneous environment. Collective flavor mixing is of interest for the
study of core-collapse supernovae, as it is thought that such oscillations may
occur at late times above the protoneutron star once the explosion has
occurred. This may have a dramatic impact on the observed neutrino signal and
nucleosynthesis in the CCSN environment [26]. Collective oscillations are
well-studied in the steady state, free streaming approximation [81]. One of
the distinct features of collective flavor mixing are the spectral swaps and
splits that arise in the neutrino and antineutrino distributions at large
radii (e.g., [82, 83]). This feature causes the neutrino spectra to deviate
from the thermal-like spectra emitted from the neutrinosphere. It is also well
established from such studies that the outcome of these calculations is
heavily dependent on fully resolving the angular distribution of neutrinos
[58]. This makes it a challenging test of a moment-based approach.

We present three simulations of this system using both the one-moment and two-
moment closures in the following subsections. The first calculation in Section
III.2.1 is a reasonably realistic representation of the conditions outside a
protoneutron star after the onset of explosion. The second calculation in
Section III.2.2 is a modified version of this same setup, where we
artificially adjust the parameters in order to probe the behavior of the
moment method when flavor transformations occur in regions where the neutrino
distribution is not highly forward-peaked. Finally, the third calculation in
Section III.2.3 is the interesting and demanding case of a double spectral
split which tests the ability of moments to track the correct degree of
coherence in the neutrinos.

a |  $L_{a}$ [ergs/s] |  $\langle q\rangle_{a}$ [MeV] |  $T_{a}$ [MeV] | $\eta_{a}$ | ${\beta_{a}}$ | $F_{r}/E$ | $P_{rr}/E$
---|---|---|---|---|---|---|---
$\nu_{e}$ | $4.1\times 10^{52}$ | 9.4 | 2.1 | 3.9 | 0.0 | 0.9975 | 0.995
$\bar{\nu}_{e}$ | $4.3\times 10^{52}$ | 13.0 | 3.5 | 2.3 | 0.0 | 0.9975 | 0.995
$\nu_{x}$ | $3.95\times 10^{51}$ | 15.8 | 4.4 | 2.1 | 0.0 | 0.9975 | 0.995
$\bar{\nu}_{x}$ | $3.95\times 10^{51}$ | 15.8 | 4.4 | 2.1 | 0.0 | 0.9975 | 0.995
Table 1: Parameters for the realistic collective oscillation simulation in
Section III.2.1. Listed are the neutrino luminosity $L_{a}$, average energy
$\langle q\rangle_{a}$, temperature $T_{a}$, and chemical potential divided by
the temperature $\eta_{a}$ used in Equation 22. These are modified from the
case of $p=10$ and $q=3.5$ (different from the $q$ in this work that
represents neutrino energy) from Table 6 in Keil  _et al._ [84], with the
electron flavor luminosities increased by a factor of 10. All angular
distributions are described by $\beta_{a}=0$ in Equation 23. The initial flux
factors $F_{r}/E$ and initial Eddington factors $P_{rr}/E$ are also provided.
Figure 4: Initial and final spectra versus energy for the test in Section
III.2.1 and Table 1. Neutrinos are in the top two panels and antineutrinos are
in the bottom two. The left column shows the two-moment simulation and the
right column shows the multi-angle simulation. The black curves are the
electron flavor and the green curves are the $x$ flavor. The solid curves show
the initial spectra and the dashed curves show the final spectra. The neutrino
spectrum in all cases shows the spectral splits common in collective mixing.
The antineutrinos show an almost complete spectral swap.

####  III.2.1 Instability far from the neutrinosphere

In Table 1 we present the values for the neutrino luminosity $L_{a}$, average
energy $\langle q\rangle_{a}$, temperature $T_{a}$, and chemical potential
$\eta_{a}$ to be used in the first collective test calculation. These values
are a modified version of the case $p=10$ and $q=3.5$ in Table 6 of Keil  _et
al._ [84], with the electron neutrino and antineutrino luminosities increased
by a factor of 10. This makes the electron neutrino and electron antineutrino
luminosities typical of what one expects during the accretion phase of an
iron-core collapse supernova - see, for example, figure 2 in [85]. The inner
computational boundary was chosen to be $100\,\mathrm{km}$ and the system was
allowed to evolve until the neutrinos reached the radius of
$400\,\mathrm{km}$. We set the mixing parameters to $\Delta
m^{2}_{12}=-2.7\times 10^{-3}\;{\rm eV}^{2}$ (similar to $\Delta m_{32}^{2}$
in the Inverted Mass ordering [86]) and $\theta_{12}=0.01$ in order to emulate
the effect of matter suppression. The neutrinosphere is set to
$R_{\nu}=10\,\mathrm{km}$ for all neutrino and antineutrino flavors.

In Figure 3 we show the results of the calculation. The figure shows The ratio
of the ‘ee’ element of luminosity matrix relative to its initial value as a
function of the radius. In the figure we see the moment calculations exhibit a
flavor instability at $126\,\mathrm{km}$, similar to, but slightly larger than
the radius of $119\,\mathrm{km}$ at which the flavor instability is seen in
the multi-angle simulation. For both approaches, once the instability has
begun we see that antineutrinos experience an almost complete flavor swap. The
figure also indicates that there is little difference between the moment
calculations that use the one-moment and two-moment closures. This similarity
was not enforced nor expected and motivated the next test problem in Section
III.2.2.

There is also remarkable agreement between the flavor-transformed spectra
between the moment and multi-angle methods. Figure 4 shows the initial and
final spectra from the one-moment, two moment, and multi-angle calculations.
In both the moment-based and multi-angle calculations, the antineutrino
spectra show an almost complete swap between the $\bar{\nu}_{e}$ and
$\bar{\nu}_{x}$ flavors and the neutrinos show a split in the spectra at about
$25\,\mathrm{MeV}$. This phenomenon has been seen many times in previous
studies starting with Duan _et al._ [87, 88].

####  III.2.2 An instability close to the neutrinosphere

a |  $L_{a}$ [ergs/s] |  $\langle q\rangle_{a}$ [MeV] |  $T_{a}$ [MeV] | $\eta_{a}$ | ${\beta_{a}}$ | $F_{r}/E$ | $P_{rr}/E$
---|---|---|---|---|---|---|---
$\nu_{e}$ | $2.050\times 10^{49}$ | 9.4 | 2.1 | 3.9 | 0.0 | 0.5 | 0.33
$\bar{\nu}_{e}$ | $2.550\times 10^{49}$ | 13.0 | 3.5 | 2.3 | 0.0 | 0.5 | 0.33
$\nu_{x}$ | $1.698\times 10^{49}$ | 15.8 | 4.4 | 2.1 | 0.0 | 0.5 | 0.33
$\bar{\nu}_{x}$ | $1.698\times 10^{49}$ | 15.8 | 4.4 | 2.1 | 0.0 | 0.5 | 0.33
Table 2: Parameters for the collective oscillation simulation discussed in
Section III.2.2 designed to expose independent evolution of multiple moments.
Listed are the luminosity $L_{a}$, average energy $\langle q\rangle_{a}$,
temperature $T_{a}$, chemical potential divided by the temperature $\eta_{a}$,
the angular distribution parameter $\beta_{a}$, the initial flux factor
$F_{r}/E$ and initial Eddington factor and $P_{rr}/E$. Due to the rapid
increase of the run time of the multi-angle calculation as the inner boundary
is moved towards the neutrinosphere, the multi-angle calculation was started
at 1 km above the neutrinosphere.

As we pointed out in Section III.2.1, the high degree of concordance between
the results using the one and two-moment closures was not enforced nor
expected, but in hindsight is perhaps not surprising. At radii well beyond the
neutrinosphere, there is very little difference between the flux and pressure
moments, since the flux and Eddington factors are close to unity. One cannot
expect to resolve fine angular features in a pencil-beam distribution using
only coarse moments. If this interpretation is correct, we might reasonably
expect to observe more significant differences between the one and two-moment
calculations when flavor instabilities occur closer to the neutrinosphere
(i.e., where the flux and Eddington factors are not close to unity). For this
reason we present results from a second self-interaction test case where we
engender flavor transformation closer to the neutrinosphere by artificially
adjusting the neutrino luminosities. The new set of parameters we use for this
calculation are shown in Table 2. In addition, for this test the inner
boundary of the computational domain is at the neutrinosphere
($R_{\nu}=10\,\mathrm{km}$).

Figure 5: The ratio of the ‘ee’ element of luminosity matrix relative to its
initial value, for the modified collective oscillation test in Section III.2.2
and Table 2. The top panel is for neutrinos and bottom panel is for
antineutrinos. The orange curve shows the one-moment simulation, blue curve
shows the two-moment simulation, and the black curve shows the multi-angle
simulation.

The evolution of the ‘ee’ element of luminosity matrix relative to its initial
value from this second self-interaction calculation are shown in Figure 5. The
system behaves similarly to the previous test problem studied in Section
III.2.1 in that the moment calculations (blue and orange) closely track each
other and exhibit an onset of flavor transformation in approximately the same
location as in the multi-angle calculations. However, due to the decreased
neutrino luminosities, the instability starts at $r\approx 14\,\mathrm{km}$,
much closer to the neutrinosphere than the onset at $r\approx
120\,\mathrm{km}$ in Section III.2.1. The moment calculations agree with the
multi-angle calculations about the onset of the instability even better than
in the previous case.

In the insets in Figure 5, we show the same ratio of the ratio of the ‘ee’
element of luminosity matrix relative to its initial value, during the first
five kilometers of the calculation. We observe that the ratios from the multi-
angle and one-moment calculations are small and slightly below unity over this
region, but that the two-moment calculation are slightly above unity until
significant flavor conversion occurs at $14\,\mathrm{km}$. Our investigation
into the origin of the greater-than-unity luminosities - which are not
physical - are presented in Appendix B. These investigations revealed that the
origin is not numerical error and are likely due to the closure.

The similarity between the results for the one- and two-moment closures shown
in Figure 5 cannot be attributed to the fact that the flux and Eddington
factors are close as in Section III.2.1. Instead, we have found that this
behavior arises because our choice of closure is so self-consistent that the
evolution equation for the flux (Equation 8) is essentially identical using
both the one-moment and two-moment closures. Thinking first about the analytic
relationships for the moments in a bulb model without flavor transformation,
one can relate the pressure to the flux as $P_{rr}=D(r)\,F_{r}$, where $D(r)$
is a scalar function given by

$\displaystyle D(r)$
$\displaystyle=\frac{2\,(1+\cos\theta_{max}+\cos^{2}\theta_{max})}{3\,(1+\cos\theta_{max})},$
(37)
$\displaystyle=\frac{2}{3\,r}\left(\frac{2\,r^{2}-R^{2}_{\nu}+r\,\sqrt{r^{2}-R^{2}_{\nu}}}{r+\sqrt{r^{2}-R^{2}_{\nu}}}\right).$

Inserting this relationship into Equations 8 and 10, the evolution equation
for the steady-state flux in the two-moment approach becomes the same as the
radial evolution of the steady state flux using the one-moment closure,
indicating the self-consistency of our choice of closures. Thus, wherever the
flux and pressure as calculated by the two-moment approach are related by
Equation 37, the evolution of the flux in the one- and two-moment systems will
evolve identically.

Figure 6: The ratio of the elements from the pressure tensor to the
corresponding elements from the flux as a function of the radius $r$ for the
$15\;{\rm MeV}$ neutrinos in the two-moment calculation using the parameters
given in Section III.2.2 and Table 2.  Figure 7: The difference between the
phase angles of the off-diagonal elements of the pressure and flux moments for
15 MeV neutrinos from the two-moment simulation (cyan) and the multi-angle
simulation (orange) in the modified collective oscillation setup of Section
III.2.2 and Table 2.

To examine whether this coherence of the flux and pressure moment actually
occurs in our calculations, in Figure 6 we plot $D(r)$ (orange) together with
the ratio of the corresponding elements from the pressure and flux moments for
the $15\;{\rm MeV}$ neutrinos. We find the ratios of
$P_{rr}^{(ee)}/F_{r}^{(ee)}$ (blue dashed) and $P_{rr}^{(xx)}/F_{r}^{(xx)}$
(green dashed, not visible under the blue dashed curve) match the function
$D(r)$ very well, indicating that the diagonal moments grow and shrink due to
flavor transformation at the same rate. The ratio
$|P_{rr}^{(ex)}|/|F_{r}^{(ex)}|$ (purple) does not follow $D(r)$ before the
instability, as the flavor off-diagonal elements are very small. Once flavor
transformation begins at $r\sim 14\;{\rm km}$, the flavor off-diagonal ratio
more closely follows $D(r)$, though not as well as the diagonal elements.

Figure 7 shows the phase angle difference
$\xi=\mathrm{arg}(P_{rr}^{(ex)})-\mathrm{arg}(F_{r}^{(ex)})$ between the off-
diagonal elements of the flux and pressure moments from the two-moment and
multi-angle calculations. The figure shows how the phase difference for the
multi-angle calculation is always small, indicating that the flux and pressure
moments experience tightly coupled evolution. The phase difference in the two-
moment calculation is initially highly variable, since the off-diagonal
components are very small. However, once the flavor transformation begins
around $r=14\;{\rm km}$, the two moments become very coherent and the phase
difference does not exceed a few degrees. Taken together, Figures 6 and 7 show
that $P_{rr}=D(r)\,F_{r}$ is a good approximation for this problem and the
coherence between the two moments reduces the system of moment equations to
that for a single moment. Note that this coherence between the pressure and
flux in the two-moment scheme is not imposed, it emerges as the calculation
proceeds.

Figure 8: Initial and final spectra versus energy for neutrinos from
simulations using parameters in Table 2 and Section III.2.2. Neutrinos are in
the top two panels and antineutrinos are in the bottom two. The left column
shows the two-moment simulation and the right column shows the multi-angle
simulation. The black lines are the electron flavor and the green lines are
the $x$ flavor. The solid lines show the initial spectra and the dashed lines
show the final spectra.

Finally, Figure 8 shows the spectra of the neutrinos and antineutrinos at
$50\;{\rm km}$ from the two-moment and multi-angle calculations. The spectra
from the two methods agrees well with only slight differences in the high-
energy tails. The moment-based approach accurately calculates the spectral
evolution as well as the gross, integrated behavior.

Figure 9: The ratio of the ‘ee’ element of luminosity matrix relative to its
initial value, for the modified collective oscillation test in Section III.2.3
and Table 3. The top panel is for neutrinos and bottom panel is for
antineutrinos. The orange curve shows the one-moment simulation, blue curve
shows the two-moment simulation, and the black curve shows the multi-angle
simulation. Figure 10: The phase of the off-diagonal element of the $H_{E}$
component of the Hamiltonian (top panel), the $H_{F}$ contribution (middle
panel) and the difference between them (bottom panel) as a function of the
radial coordinate. Blue lines are the results from the moment-based approach
using the two-moment closure, the black are from the multi-angle calculation.
Figure 11: Initial and final spectra versus energy for neutrinos from
simulations using parameters in Table 3 and Section III.2.3. Neutrinos are in
the top two panels and antineutrinos are in the bottom two. The left column
shows the two-moment simulation and the right column shows the multi-angle
simulation. The black lines are the electron flavor and the green lines are
the $x$ flavor. The solid lines show the initial spectra and the dashed lines
show the final spectra.

####  III.2.3 A case with multiple spectral splits

a |  $L_{a}$ [ergs/s] |  $\langle q\rangle_{a}$ [MeV] |  $T_{a}$ [MeV] | $\eta_{a}$ | ${\beta_{a}}$ | $F_{r}/E$ | $P_{rr}/E$
---|---|---|---|---|---|---|---
$\nu_{e}$ | $1.8\times 10^{52}$ | 12.0 | 2.1 | 3.9 | 0.0 | 0.933 | 0.872
$\bar{\nu}_{e}$ | $2.2\times 10^{52}$ | 15.0 | 3.5 | 2.3 | 0.0 | 0.933 | 0.872
$\nu_{x}$ | $2.7\times 10^{52}$ | 18.0 | 4.4 | 2.1 | 0.0 | 0.933 | 0.872
$\bar{\nu}_{x}$ | $2.7\times 10^{52}$ | 18.0 | 4.4 | 2.1 | 0.0 | 0.933 | 0.872
Table 3: Parameters for the collective oscillation simulation discussed in
Section III.2.3 which lead to a case of multiple spectral splits. Listed are
the luminosity $L_{a}$, average energy $\langle q\rangle_{a}$, temperature
$T_{a}$, chemical potential divided by the temperature $\eta_{a}$, the angular
distribution parameter $\beta_{a}$, the initial flux factor $F_{r}/E$ and
initial Eddington factor and $P_{rr}/E$.

Our final test of the moment code is for a case which leads to multiple
spectral splits. The parameters for this test problem as shown in table 3 and
are a slight adjustment from those used in the first example in [88]. The
luminosities and mean energies are such that the number fluxes of all four
flavors are almost equal, with the electron antineutrino number flux $\sim
2\%$ smaller than the others. This near-equality of the number fluxes leads to
a small self-interaction potential which tends to move the instability close
to the proto-neutron star.

Figure 9 shows the evolution of the ‘ee’ element of the neutrino and
antineutrino luminosity matrices as a function of the radius. A flavor
instability is seen to occur around $r\sim 40\;{\rm km}$ and, once again, both
moment calculations with the two different closures and the multi-angle
calculation are in good agreement about the location of the instability. For
the next $\sim 20\;{\rm km}$ thereafter the three calculations track one
another closely but at $r\sim 60\;{\rm km}$ there is a noticeable change in
the multi-angle calculation compared to the moment-based approach. Beyond
$r\sim 60\;{\rm km}$ the multi-angle results exhibit rapid oscillations which
decrease in amplitude so that by $r\sim 200\;{\rm km}$ they are no longer
observable. The moment-based approaches change noticeably at $r\sim 100\;{\rm
km}$ and while the oscillations in the elements of the luminosity matrix also
decrease with distance, it is over a much longer scale. In contrast to the two
previous self-interaction test problems, there are substantial difference of
the luminosities at large radii of the different approaches.

We have explored why the multi-angle and moment based results differ in this
test problem by focusing upon the coherence of the two contributions to the
self-interaction. Figure 10 shows the phase angles $\xi_{E}$ and $\xi_{F}$ of
the off-diagonal elements of the two contributions to the self-interaction,
$H_{E}$ and $H_{F}$, respectively (i.e. $\xi_{E}=\arg(H_{E}^{e\mu})$,
$\xi_{E}=\arg(H_{F}^{e\mu})$ ) and the difference between them for the moment-
based approach using the two-moment closure, and from the multi-angle
calculation. The evolution of these phase angles separately with radius
indicate the off-diagonal elements of $H_{E}$ and $H_{F}$ rotate rapidly in
the Argand plane, but even so, for $r\lesssim 60\;{\rm km}$ the difference
between the phase angles is miniscule and the two terms are coherent. Beyond
$r\sim 60\;{\rm km}$ the evolution of the phase angles for the multi-angle
calculation changes noticeably. At first glance the change in behavior of each
phase angle appears to be similar but closer inspection reveals this not to be
the case and the difference between the two phases becomes non-zero indicating
that actually the degree of coherence between the two contributions to
$H_{SI}$ is weaker. The loss of coherence of these two contributions to the
off-diagonal elements of the self-interaction Hamiltonian reduces the amount
of flavor transformation in the multi-angle calculations beyond $r\sim
60\;{\rm km}$. Fluctuations in $\xi_{F}-\xi_{E}$ are not seen in the moment-
based results and thus the two contributions to $H_{SI}$ remain coherent
leading to significant flavor transformation after $r\sim 60\;{\rm km}$. In
order to verify this difference between the two approaches we have repeated
the calculations using different numbers of energy and angle bins, and even a
different code, and found the same results every time. Thus we suspect the
loss of coherence in the multi-angle calculation is due to the significant
cancellation in the self-interaction Hamiltonian due to the spectral
parameters chosen for this test problem. Less cancellation - as in the
previous two self-interaction test problems - leads to stronger coherence and
better agreement between the moment-based and multi-angle results.

The differences between the two approaches is also seen in the spectra shown
in Figure 11. The spectra from the moment calculation shows two splits in the
neutrinos at $E\sim 4\;{\rm MeV}$ and $E\sim 30\;{\rm MeV}$, and $E\sim
25\;{\rm MeV}$ in the antineutrinos, with complete swaps of the spectra
outside the split regions. In contrast, the spectra from the multi-angle
calculation are very close to a 50:50 mixture of the initial spectra albeit
with observable changes in the exact amount of mixing at the same split
energies seen more clearly in the moment-based approach. These results are a
close match to those seen in [88] where single-angle calculations were seen to
produce sharp splits and complete swaps, while the splits in the multi-angle
calculations were less sharp and the swaps incomplete.

##  IV Discussion & Conclusions

We have presented here a method for modeling neutrino flavor mixing using one
and two moments. We then investigated how well it captures the mixing effects
of analytic predictions and of more exact but more computationally expensive
multi-angle calculations by considering several problems in a bulb-model
geometry representative of core-collapse supernovae. Overall, the moment-based
method reproduces results from a multi-angle method surprisingly well.
However, errors can emerge in certain circumstances as a result of the
tendency of the moment method to maintain an artificially high level of
coherence between the simulated moments when using a scalar closure
independent of the level at which we truncate the evolution equation tower.
This supports the claim that the largest source of error is the nature of the
closure and not the number of evolved moments.

The first problem was the flavor evolution of neutrinos emitted from a
neutrino bulb in matter with a density chosen to put the 1 MeV neutrinos on
the MSW resonance for the given mixing parameters. We found that the moment
calculations overestimated the amplitude of the flavor mixing at a given
radial point because they were not able to accurately account for the
incoherence of neutrinos which had traveled along different trajectories. As
the emission at the neutrinosphere becomes more forward-peaked i.e. a greater
proportion of the neutrinos are emitted along trajectories close to the radial
direction, the transition probabilities from the multi-angle calculation
become more similar with those from the moment calculation.

The last three test problems were used to study the case of collective flavor
transformation due to neutrino self interaction. In the first of these
calculations, designed with a reasonably realistic inner boundary condition,
the moment and multi-angle results were in good agreement in that the onset of
the transformation differed by just a few kilometers. The results from the
one-moment and two-moment calculations were almost identical. When we adjusted
the parameters so as to produce flavor transformation much closer to the
neutrinosphere (i.e., where the flux and Eddington factors are much lower), we
again found the multi-angle and moment calculations were in good agreement on
where the flavor transformation begins and that, once again, the one-moment
and two-moment calculations yield essentially the same radial evolution of the
luminosities. We then demonstrated that the similarity of the two moment-based
approaches is due to a sustained synchronization between the flux and pressure
moments which effectively collapses the tower of moment equations to that for
a single moment. This coherence between the flux and pressure moments appears
to emerge naturally in steady state self-interaction situations. As a final
attempt to elicit different outcomes from different methods, we simulate
conditions constructed to produce multiple spectral swaps. Once again, the
moment method correctly predicted the onset of instability and for some range
thereafter, the two approaches were in strong agreement. But eventually the
two approaches diverged due to the artificially high levels of coherence in
the moment approach compared to the multi-angle and this, in this case, the
difference resulted in much larger errors.

Overall, the various tests we undertook indicate that a moment-based approach
does well at capturing the overall neutrino flavor transformation seen in the
more computationally expensive, multi-angle, calculation. The onset of flavor
mixing is predicted with an accuracy of a few kilometers and the spectra are
similar. However, we do find differences between the two approaches indicating
further studies are required before we can completely enjoy the benefits of
the moment-based approach. Our second self-interaction test problem revealed
that the moment-based approaches can yield, at least temporarily, unphysical
solutions, and we also saw in our third self-interaction test case how over-
estimation of the coherence can lead to different final spectra. The
differences between the multi-angle and moment-based results are entirely due
to the closure. The closures we have used in this paper are scalar relations,
but for quantum moments in general, one would expect phase differences between
the off-diagonal elements of the two moments used in the closure. A scalar
closure like those used here cannot generate such a phase difference. To
permit the phase differences to appear, one requires a more general, quantum
closure. We shall explore quantum closures in future work.

In summary, our results indicate that moment-based schemes are an inexpensive
approach to neutrino transport that are able to capture most of the flavor
mixing phenomenology seen in more-exact, but more expensive, calculations. We
caution that any enthusiasm for a moment-based approach to neutrino
oscillations must be tempered by remembering that that any approach which
truncates the tower of moment equations cannot converge to those from the full
quantum kinetic equations and thus moment-based approaches will continue to
require verification against those derived from less approximate methods.
Should further, more demanding comparisons of moments and less approximate
methods, (especially comparisons which do not assume a steady state) reveal
that moment-based methods perform well in those situations, too, we would
proffer moment methods are a promising and viable avenue for including
neutrino transformation in hydrodynamical simulations.

###### Acknowledgements.

The authors would like to thank S. Couch and E. P. O’Connor for their advice
and support through this project. MLW was supported by a NSF Astronomy and
Astrophysics Postdoctoral Fellowship under award AST-1801844. TC, MM, JPK,
GCM, EG, and CF are supported by the DOE Office of Nuclear Physics award DE-
FG02-02ER41216. E. G. acknowledges support by the National Science Foundation
grant No. PHY-1430152 (JINA Center for the Evolution of the Elements). This
work was partially enabled by the National Science Foundation under grants No.
PHY-1630782 (N3AS hub) and No. PHY-2020275 (N3AS PFC). SR was supported by a
NSF Astronomy and Astrophysics Postdoctoral Fellowship under award
AST-2001760. The calculations presented in this paper were undertaken with the
support of Michigan State University through computational resources provided
by the Institute for Cyber-Enabled Research, and on the _Payne_ machine at
NCSU which is supported in part by the Research Corporation for Science
Advancement. The authors would like to acknowledge the use of the following
software: Matplotlib [89], Numpy [90], and SciPy [91].

Figure 12: The difference in the transition probability as a function of
radius for the 15 MeV neutrinos from a two-moment calculation with an energy
resolution of $50\;{\rm keV}$ compared to the fiducial energy resolution of
$100\;{\rm keV}$. The calculation uses the parameters shown in Table 2.

##  Appendix A SQA

For the multi-angle calculations we use SQA, which is an implementation of the
Bulb Model described in Duan  _et al._ [92]. The momenta of the neutrinos and
antineutrinos at the neutrinosphere are discretized into $N_{E}$ energy bins
and $N_{A}$ angle bins. The energy resolution is uniform, but for the
direction we adopt uniform resolution in the quantity $u=\sin^{2}\theta_{R}$
where $\theta_{R}$ is the angle at which the neutrino was emitted at the
neutrinosphere relative to the radial direction. This distribution gives
greater weight to the angles emitted at large angles $\theta_{R}$, since these
rays have the largest dispersion in path lengths from the neutrinosphere to a
given radial point. We introduce an evolution matrix $S$ that relates the
initial density matrix for neutrinos moving with each momentum $\vec{p}_{0}$
at the neutrinosphere to the density matrix for neutrinos moving at some
distance $\lambda$ along the same trajectory with momentum $\vec{p}(\lambda)$.
That is,

$\rho(\lambda,\vec{p}_{0})=S(\lambda,\vec{p}_{0})\,\rho(R_{\nu},\vec{p}_{0})\,S^{\dagger}(\lambda,\vec{p}_{0}).$
(38)

A matrix $\bar{S}(\lambda,\vec{p}_{0})$ plays the same role for antineutrinos.
The matrix $S(\lambda,\vec{p}_{0})$ evolves according to the Schrödinger
equation

$i\frac{dS}{d\lambda}=H\left(\lambda,\vec{p}\right)S.$ (39)

For the antineutrinos, the matrix $\bar{S}$ similarly evolves according to the
Hamiltonian $\bar{H}$. As in the moment method described in Section II, the
Hamilton is composed of three parts: vacuum $H_{V}$, matter $H_{M}$, and self-
interaction $H_{\mathrm{SI}}$. The vacuum and matter terms were given
previously in Equations 14 and 15, and the self interaction term is

$H_{\mathrm{SI}}\left(r,{\vec{p}}\right)=\sqrt{2}\,G_{F}\int{\left({1-{\bf{\hat{p}}}\cdot{\bf{\hat{q}}}}\right)\left[{\rho(r,\vec{q})-\bar{\rho}^{*}(r,\vec{q})}\right]q^{2}dq\,d\Omega}.$
(40)

The density matrices are evolved to discrete radial points $r(\lambda)$ common
among all trajectories so the Hamiltonian can be straightforwardly integrated.

To ensure unitarity, the matrices $S$ and $\bar{S}$ are parameterized by four
real variables: three are the angles defining a unit vector in a four-
dimensional Euclidian flavor space, and the fourth is the phase of the
determinant. The differential equations for all of the parameters are solved
simultaneously with an explicit Runge-Kutta integrator that uses an adaptive
step size and the Cash-Karp parameter set. Several sources of numerical errors
in multi-angle codes have been identified over the years [93, 57, 58, 94]. The
numerical accuracy of the SQA calculations presented in this paper is
estimated to be less than 1% based on the lack of visible differences between
the results using different numbers of energy or angle bins, and initial
radii. A related code, IsotropicSQA, was used for solving the neutrino quantum
kinetic equations in isotropic and homogeneous conditions [95] is available at
[71].

##  Appendix B Numerical Errors in the Moment-Based Calculations

During our examination of the results from the moment-based code using the
two-moment closure, we observed that the luminosites became slightly and
temporarily greater than their initial values, as shown in the inset of Figure
5. The same ratio from the code using the one-moment closure scheme did not
contain this feature. We have explored the reason for the appearance of these
unphysical results from the two-moment scheme and attempted several strategies
to remove them. In this appendix we describe those efforts.

The greater-than-initial luminosities are present only when the self-
interaction contribution to the Hamiltonian is included, and they occur for
both neutrinos and antineutrinos. We have no reason to suspect that the errors
are numerical in origin generated from the ODE integrator we used, since
decreasing the error tolerance for the integrator does not improve the
solution. We also changed the algorithm from an explicit second-order
‘midpoint’ Runge-Kutta integrator to an adaptive step size routine using the
Cash-Karp parameter set and a fractional error tolerance per step of
$10^{-10}$, and found the negative factors persisted. We also tried changing
the parameterization of the moment matrices by writing the diagonal elements
as $M^{(ee)}=Tr(M)\,\cos^{2}\theta_{M}$ and
$M^{(xx)}=Tr(M)\,\sin^{2}\theta_{M}$ where $Tr(M)$ is the trace and
$\theta_{M}$ an angle, in an effort to prohibit the unphysical behavior of the
diagonal elements. However, when we ran our code using the two-moment closure
to evolve these new parameters, the code terminated prematurely when the step
size shrank to zero at the point where the original version crossed into the
unphysical solution space.

Figure 13: The transition probability as a function of radius as computed from
a single-energy calculation using the parameters shown in Table 2. The energy
of the neutrinos and antineutrino is chosen to be $5\;{\rm MeV}$ in order to
produce an instability in approximately the same location seen in the multi-
energy calculation.

The greater-than-initial luminosities are not dependent upon the number of
energy bins used in the calculation. While the results in Section III.2.2 use
and energy resolution of $100\;{\rm keV}$, we repeated the calculation using a
energy resolution of $50\;{\rm keV}$. The evolution of the transition
probability (the transition probability is defined in Eq. 34) for the
$15\;{\rm MeV}$ neutrinos from the two calculations is shown in Figure 12. We
observe changes in the transition probability that are of order 0.001% in the
interval $10\;{\rm km}\leq r\leq 14\;{\rm km}$ where the greater-than-initial
luminosities appear. We also ran a single energy calculation - for which
energy resolution is not a factor - using an energy of $5\;{\rm MeV}$. The
choice of $5\;{\rm MeV}$ produces an onset of flavor transformation at
$r\approx 13.5\,\mathrm{km}$, similar to that seen in the multi-energy
calculation. The results from this single-energy calculation are shown in
Figure (13) and we observe negative transition probabilities. Thus we do not
believe numerical error from the energy resolution to be the origin of the
negative factors.

We suspect the origin of the errors to be a combination of using moments, and
the closure. In general moments do not evolve according to a unitary operator
because a moment is an integrated quantity which aggregates information about
the neutrinos traveling along separate trajectories. This is reflected in the
evolution equations themselves. The last term on the right hand side of the
moment evolution equations (Equations 8-11) unitarily evolves the moment’s
flavor structure and cannot lead to negative probabilities except for
numerical error. The geometrical term (the last term on the left hand side of
the equation) leads to non-unitary evolution but this evolution is simply a
scaling and cannot lead to the problems we see. It is the first commutator on
the right hand side which is the source of physical non-unitary evolution of
the moments and depends on the choice of closure. Thus we suspect that the
closure we use, while geometrically justified and always realizable, may not
be sufficient to enforce physical evolution of the moment equations. An
analysis of various options for closures which enforce that the solution
always remain physical is beyond the scope of this paper.

## References

  * Colgate and White [1966] S. A. Colgate and R. H. White, The Hydrodynamic Behavior of Supernovae Explosions, The Astrophysical Journal 143, 626 (1966).
  * Bethe [1990] H. A. Bethe, Supernova mechanisms, Reviews of Modern Physics 62, 801 (1990), publisher: American Physical Society.
  * Kotake  _et al._ [2006] K. Kotake, K. Sato, and K. Takahashi, Explosion mechanism, neutrino burst and gravitational wave in core-collapse supernovae, Reports on Progress in Physics 69, 971 (2006).
  * Janka  _et al._ [2007] H. T. Janka, K. Langanke, A. Marek, G. Martínez-Pinedo, and B. Müller, Theory of core-collapse supernovae, Physics Reports The Hans Bethe Centennial Volume 1906-2006, 442, 38 (2007).
  * Janka [2012] H.-T. Janka, Explosion Mechanisms of Core-Collapse Supernovae, Annual Review of Nuclear and Particle Science 62, 407 (2012).
  * Burrows [2013] A. Burrows, Colloquium: Perspectives on core-collapse supernova theory, Reviews of Modern Physics 85, 245 (2013), publisher: American Physical Society.
  * Foglizzo  _et al._ [2015] T. Foglizzo, R. Kazeroni, J. Guilet, F. Masset, M. González, B. K. Krueger, J. Novak, M. Oertel, J. Margueron, J. Faure, N. Martin, P. Blottiau, B. Peres, and G. Durand, The Explosion Mechanism of Core-Collapse Supernovae: Progress in Supernova Theory and Experiments, Publications of the Astronomical Society of Australia 32, 10.1017/pasa.2015.9 (2015), publisher: Cambridge University Press.
  * Lentz  _et al._ [2015] E. J. Lentz, S. W. Bruenn, W. R. Hix, A. Mezzacappa, O. E. B. Messer, E. Endeve, J. M. Blondin, J. A. Harris, P. Marronetti, and K. N. Yakunin, THREE-DIMENSIONAL CORE-COLLAPSE SUPERNOVA SIMULATED USING A 15 M $ødot$ PROGENITOR, The Astrophysical Journal 807, L31 (2015).
  * Müller  _et al._ [2017] B. Müller, T. Melson, A. Heger, and H.-T. Janka, Supernova simulations from a 3D progenitor model – Impact of perturbations and evolution of explosion properties, Monthly Notices of the Royal Astronomical Society 472, 491 (2017).
  * Ott  _et al._ [2018] C. D. Ott, L. F. Roberts, A. d. S. Schneider, J. M. Fedrow, R. Haas, and E. Schnetter, The Progenitor Dependence of Core-collapse Supernovae from Three-dimensional Simulations with Progenitor Models of 12–40 M $ødot$, The Astrophysical Journal 855, L3 (2018).
  * Summa  _et al._ [2018] A. Summa, H.-T. Janka, T. Melson, and A. Marek, Rotation-supported Neutrino-driven Supernova Explosions in Three Dimensions and the Critical Luminosity Condition, The Astrophysical Journal 852, 28 (2018), publisher: American Astronomical Society.
  * O’Connor and Couch [2018] E. P. O’Connor and S. M. Couch, Exploring Fundamentally Three-dimensional Phenomena in High-fidelity Simulations of Core-collapse Supernovae, The Astrophysical Journal 865, 81 (2018).
  * Nakamura  _et al._ [2019] K. Nakamura, T. Takiwaki, and K. Kotake, Long-term simulations of multi-dimensional core-collapse supernovae: Implications for neutron star kicks, Publications of the Astronomical Society of Japan 71, 10.1093/pasj/psz080 (2019).
  * Glas  _et al._ [2019] R. Glas, H. T. Janka, T. Melson, G. Stockinger, and O. Just, Effects of LESA in Three-dimensional Supernova Simulations with Multidimensional and Ray-by-ray-plus Neutrino Transport, The Astrophysical Journal 881, 36 (2019).
  * Glas [2019] R. Glas, Three-dimensional Core-collapse Supernova Simulations with Multidimensional Neutrino Transport Compared to the Ray-by-ray-plus Approximation, The Astrophysical Journal 873, 10.3847/1538-4357/ab0423 (2019).
  * Burrows  _et al._ [2020] A. Burrows, D. Radice, D. Vartanyan, H. Nagakura, M. A. Skinner, and J. C. Dolence, The overarching framework of core-collapse supernova explosions as revealed by 3D fornax simulations, Monthly Notices of the Royal Astronomical Society 491, 2715 (2020).
  * Yasin  _et al._ [2020] H. Yasin, S. Schäfer, A. Arcones, and A. Schwenk, Equation of State Effects in Core-Collapse Supernovae, Phys. Rev. Lett. 124, 092701 (2020).
  * da Silva Schneider  _et al._ [2020] A. da Silva Schneider, E. O’Connor, E. Granqvist, A. Betranhandy, and S. M. Couch, Equation of State and Progenitor Dependence of Stellar-mass Black Hole Formation, The Astrophysical Journal 894, 4 (2020), arXiv:2001.10434 [astro-ph.HE] .
  * Ghosh  _et al._ [2021] S. Ghosh, N. Wolfe, and C. Fröhlich, PUSHing core-collapse supernovae to explosions in spherical symmetry V: Equation of state dependency of explosion properties, nucleosynthesis yields, and compact remnants, arXiv e-prints , arXiv:2107.13016 (2021), arXiv:2107.13016 [astro-ph.HE] .
  * Radice  _et al._ [2016] D. Radice, C. D. Ott, E. Abdikamalov, S. M. Couch, R. Haas, and E. Schnetter, Neutrino-driven convection in core-collapse supernovae: high-resolution simulations, The Astrophysical Journal 820, 76 (2016).
  * Richers  _et al._ [2017] S. Richers, H. Nagakura, C. D. Ott, J. Dolence, K. Sumiyoshi, and S. Yamada, A Detailed Comparison of Multidimensional Boltzmann Neutrino Transport Methods in Core-collapse Supernovae, The Astrophysical Journal 847, 133 (2017).
  * O’Connor  _et al._ [2018] E. O’Connor, R. Bollig, A. Burrows, S. Couch, T. Fischer, H.-T. Janka, K. Kotake, E. J. Lentz, M. Liebendörfer, O. E. B. Messer, A. Mezzacappa, T. Takiwaki, and D. Vartanyan, Global comparison of core-collapse supernova simulations in spherical symmetry, Journal of Physics G Nuclear Physics 45, 104001 (2018).
  * Just  _et al._ [2018] O. Just, R. Bollig, H.-T. Janka, M. Obergaulinger, R. Glas, and S. Nagataki, Core-collapse supernova simulations in one and two dimensions: comparison of codes and approximations, Monthly Notices of the Royal Astronomical Society 481, 4786 (2018).
  * Dasgupta  _et al._ [2012] B. Dasgupta, E. P. O’Connor, and C. D. Ott, Role of collective neutrino flavor oscillations in core-collapse supernova shock revival, Physical Review D 85, 065008 (2012).
  * Stapleford  _et al._ [2020] C. J. Stapleford, C. Fröhlich, and J. P. Kneller, Coupling neutrino oscillations and simulations of core-collapse supernovae, Physical Review D 102, 081301 (2020).
  * Duan  _et al._ [2010] H. Duan, G. M. Fuller, and Y.-Z. Qian, Collective Neutrino Oscillations, Annual Review of Nuclear and Particle Science 60, 569 (2010).
  * Bellini  _et al._ [2014] G. Bellini, L. Ludhova, G. Ranucci, and F. L. Villante, Neutrino Oscillations, Advances in High Energy Physics 2014, 1 (2014).
  * Mikheyev and Smirnov [1985] S. P. Mikheyev and A. Y. Smirnov, Resonance enhancement of oscillations in matter and solar neutrino spectroscopy, Yad. Fiz. 42, 1441 (1985), (Sov. J. Nucl. Phys. 42 913).
  * Mikheyev and Smirnov [1986] S. P. Mikheyev and A. Y. Smirnov, Neutrino oscillations in matter with varying density, in  _’86 Massive Neutrinos in Astrophysics and in Particle Physics_ , edited by O. Frackler and J. Trân Thanh Vân (Editions Frontières, Gif-sur-Yvette, 1986) p. 355.
  * Wolfenstein [1978] L. Wolfenstein, Neutrino oscillations in matter, Phys. Rev. D 17, 2369 (1978).
  * Kneller  _et al._ [2008] J. P. Kneller, G. C. McLaughlin, and J. Brockman, Oscillation Effects and Time Variation of the Supernova Neutrino Signal, Phys. Rev. D 77, 045023 (2008), arXiv:0705.3835 [astro-ph] .
  * Gava  _et al._ [2009] J. Gava, J. Kneller, C. Volpe, and G. C. McLaughlin, A Dynamical collective calculation of supernova neutrino signals, Phys. Rev. Lett. 103, 071101 (2009), arXiv:0902.0317 [hep-ph] .
  * Patton  _et al._ [2014] K. M. Patton, J. P. Kneller, and G. C. McLaughlin, Stimulated Neutrino Transformation Through Turbulence, Phys. Rev. D 89, 073022 (2014), arXiv:1310.5643 [hep-ph] .
  * Chakraborty  _et al._ [2011] S. Chakraborty, T. Fischer, A. Mirizzi, N. Saviano, and R. Tomàs, No Collective Neutrino Flavor Conversions during the Supernova Accretion Phase, Physical Review Letters 107, 151101 (2011), publisher: American Physical Society.
  * Duan and Friedland [2011] H. Duan and A. Friedland, Self-Induced Suppression of Collective Neutrino Oscillations in a Supernova, Physical Review Letters 106, 091101 (2011).
  * Stapleford  _et al._ [2016] C. J. Stapleford, D. J. Väänänen, J. P. Kneller, G. C. McLaughlin, and B. T. Shapiro, Nonstandard Neutrino Interactions in Supernovae, Phys. Rev. D 94, 093007 (2016), arXiv:1605.04903 [hep-ph] .
  * Sawyer [2005] R. F. Sawyer, Speed-up of neutrino transformations in a supernova environment, Physical Review D 72, 045003 (2005).
  * Sawyer [2016] R. F. Sawyer, Neutrino cloud instabilities just above the neutrino sphere of a supernova, Physical Review Letters 116, 081101 (2016).
  * Chakraborty  _et al._ [2016] S. Chakraborty, R. S. Hansen, I. Izaguirre, and G. Raffelt, Self-induced neutrino flavor conversion without flavor mixing, Journal of Cosmology and Astroparticle Physics 2016 (03), 042.
  * Izaguirre  _et al._ [2017] I. Izaguirre, G. Raffelt, and I. Tamborra, Fast Pairwise Conversion of Supernova Neutrinos: A Dispersion Relation Approach, Physical Review Letters 118, 021101 (2017).
  * Capozzi  _et al._ [2017] F. Capozzi, B. Dasgupta, E. Lisi, A. Marrone, and A. Mirizzi, Fast flavor conversions of supernova neutrinos: Classifying instabilities via dispersion relations, Physical Review D 96, 043016 (2017).
  * Dasgupta and Sen [2018] B. Dasgupta and M. Sen, Fast neutrino flavor conversion as oscillations in a quartic potential, Physical Review D 97, 023017 (2018).
  * Abbar and Duan [2018] S. Abbar and H. Duan, Fast neutrino flavor conversion: Roles of dense matter and spectrum crossing, Physical Review D 98, 043014 (2018).
  * Johns  _et al._ [2020a] L. Johns, H. Nagakura, G. M. Fuller, and A. Burrows, Neutrino oscillations in supernovae: angular moments and fast instabilities, Phys. Rev. D 101, 043009 (2020a), arXiv:1910.05682 [hep-ph] .
  * Capozzi  _et al._ [2021] F. Capozzi, S. Abbar, R. Bollig, and H. T. Janka, Fast neutrino flavor conversions in one-dimensional core-collapse supernova models with and without muon creation, Phys. Rev. D 103, 063013 (2021), arXiv:2012.08525 [astro-ph.HE] .
  * Johns and Nagakura [2021] L. Johns and H. Nagakura, Fast flavor instabilities and the search for neutrino angular crossings, Phys. Rev. D 103, 123012 (2021), arXiv:2104.04106 [hep-ph] .
  * Nagakura  _et al._ [2021] H. Nagakura, L. Johns, A. Burrows, and G. M. Fuller, Where, when and why: occurrence of fast-pairwise collective neutrino oscillation in three-dimensional core-collapse supernova models, Phys.Rev.D 104, 083025 (2021), arXiv:2108.07281 [astro-ph] .
  * Nagakura and Johns [2021] H. Nagakura and L. Johns, Constructing angular distributions of neutrinos in core-collapse supernovae from zeroth and first moments calibrated by full Boltzmann neutrino transport, Phys. Rev. D 103, 123025 (2021), arXiv:2104.05729 [astro-ph.HE] .
  * Abbar and Volpe [2019] S. Abbar and M. C. Volpe, On fast neutrino flavor conversion modes in the nonlinear regime, Physics Letters B 790, 545 (2019).
  * Shalgar and Tamborra [2019] S. Shalgar and I. Tamborra, On the Occurrence of Crossings between the Angular Distributions of Electron Neutrinos and Antineutrinos in the Supernova Core, The Astrophysical Journal 883, 80 (2019).
  * M. Delfan Azari  _et al._ [2020] M. Delfan Azari, S. Yamada, T. Morinaga, H. Nagakura, S. Furusawa, A. Harada, H. Okawa, W. Iwakami, and K. Sumiyoshi, Fast collective neutrino oscillations inside the neutrino sphere in core-collapse supernovae, Physical Review D 101, 023018 (2020).
  * Glas  _et al._ [2020] R. Glas, H. Thomas Janka, F. Capozzi, M. Sen, B. Dasgupta, A. Mirizzi, and G. Sigl, Fast neutrino flavor instability in the neutron-star convection layer of three-dimensional supernova models, Physical Review D 101, 063001 (2020).
  * Morinaga  _et al._ [2020] T. Morinaga, H. Nagakura, C. Kato, and S. Yamada, Fast neutrino-flavor conversion in the preshock region of core-collapse supernovae, Physical Review Research 2, 012046(R) (2020).
  * Mangano  _et al._ [2014] G. Mangano, A. Mirizzi, and N. Saviano, Damping the neutrino flavor pendulum by breaking homogeneity, Phys. Rev. D 89, 073017 (2014), arXiv:1403.1892 [hep-ph] .
  * Abbar [2021] S. Abbar, Turbulence Fingerprint on Collective Oscillations of Supernova Neutrinos, Phys. Rev. D 103, 045014 (2021), arXiv:2007.13655 [astro-ph.HE] .
  * Sigl [2021] G. Sigl, Simulations of Fast Neutrino Flavor Conversions with Interactions in Inhomogeneous Media, arXiv e-prints , arXiv:2109.00091 (2021), arXiv:2109.00091 [hep-ph] .
  * Duan  _et al._ [2008] H. Duan, G. M. Fuller, and J. Carlson, Simulating nonlinear neutrino flavor evolution, Computational Science and Discovery 1, 015007 (2008), arXiv:0803.3650 [astro-ph] .
  * Sarikas  _et al._ [2012] S. Sarikas, D. de Sousa Seixas, and G. Raffelt, Spurious instabilities in multiangle simulations of collective flavor conversion, Phys.Rev.D 86, 125020 (2012), arXiv:1210.4557 [hep-ph] .
  * Iwakami  _et al._ [2020] W. Iwakami, H. Okawa, H. Nagakura, A. Harada, S. Furusawa, K. Sumiyoshi, H. Matsufuru, and S. Yamada, Simulations of the Early Postbounce Phase of Core-collapse Supernovae in Three-dimensional Space with Full Boltzmann Neutrino Transport, The Astrophysical Journal 903, 82 (2020), publisher: American Astronomical Society.
  * Jiang  _et al._ [2014] Y.-F. Jiang, J. M. Stone, and S. W. Davis, AN ALGORITHM FOR RADIATION MAGNETOHYDRODYNAMICS BASED ON SOLVING THE TIME-DEPENDENT TRANSFER EQUATION,  213, 7 (2014), publisher: American Astronomical Society.
  * Cimerman  _et al._ [2017] N. P. Cimerman, R. Kuiper, and C. W. Ormel, Hydrodynamics of embedded planets’ first atmospheres – III. The role of radiation transport for super-Earth planets, Monthly Notices of the Royal Astronomical Society 471, 4662 (2017).
  * Birindelli  _et al._ [2017] G. Birindelli, J.-L. Feugeas, J. Caron, B. Dubroca, G. Kantor, J. Page, T. Pichard, V. Tikhonchuk, and P. Nicolaï, High performance modelling of the transport of energetic particles for photon radiotherapy, Physica Medica 42, 305 (2017).
  * Hopkins and Grudić [2019] P. F. Hopkins and M. Y. Grudić, Numerical problems in coupling photon momentum (radiation pressure) to gas, Monthly Notices of the Royal Astronomical Society 483, 4187 (2019).
  * Bloch  _et al._ [2021] H. Bloch, P. Tremblin, M. González, T. Padioleau, and E. Audit, A high-performance and portable asymptotic preserving radiation hydrodynamics code with the M1 model, Astronomy & Astrophysics 646, A123 (2021), publisher: EDP Sciences.
  * Fuksman  _et al._ [2021] J. D. M. Fuksman, H. Klahr, M. Flock, and A. Mignone, A Two-moment Radiation Hydrodynamics Scheme Applicable to Simulations of Planet Formation in Circumstellar Disks,  906, 78 (2021), publisher: American Astronomical Society.
  * Vlasenko  _et al._ [2014] A. Vlasenko, G. M. Fuller, and V. Cirigliano, Neutrino quantum kinetics, Physical Review D 89, 105004 (2014).
  * Volpe [2015] C. Volpe, Neutrino quantum kinetic equations, International Journal of Modern Physics E 24, 1541009 (2015).
  * Blaschke and Cirigliano [2016] D. N. Blaschke and V. Cirigliano, Neutrino quantum kinetic equations: The collision term, Physical Review D 94, 033009 (2016).
  * Strack and Burrows [2005] P. Strack and A. Burrows, Generalized Boltzmann formalism for oscillating neutrinos, Physical Review D 71, 093004 (2005), publisher: American Physical Society.
  * Duan and Shalgar [2014] H. Duan and S. Shalgar, Multipole expansion method for supernova neutrino oscillations, JCAP 2014 (10), 084, arXiv:1407.7861 [hep-ph] .
  * Richers  _et al._ [2019] S. A. Richers, G. C. McLaughlin, J. P. Kneller, and A. Vlasenko, Neutrino quantum kinetics in compact objects, Physical Review D 99, 123014 (2019).
  * Johns  _et al._ [2020b] L. Johns, H. Nagakura, G. M. Fuller, and A. Burrows, Fast oscillations, collisionless relaxation, and spurious evolution of supernova neutrino flavor, Physical Review D 102, 103017 (2020b).
  * Zhang and Burrows [2013] Y. Zhang and A. Burrows, Transport equations for oscillating neutrinos, Physical Review D 88, 105009 (2013).
  * Thorne [1981] K. S. Thorne, Relativistic radiative transfer - Moment formalisms, Monthly Notices of the Royal Astronomical Society 194, 439 (1981).
  * Smit  _et al._ [2000] J. M. Smit, L. J. van den Horn, and S. A. Bludman, Closure in flux-limited neutrino diffusion and two-moment transport, Astronomy and Astrophysics 356, 559 (2000).
  * Murchikova  _et al._ [2017] E. M. Murchikova, E. Abdikamalov, and T. Urbatsch, Analytic closures for M1 neutrino transport, Monthly Notices of the Royal Astronomical Society 469, 1725 (2017).
  * Hayes and Norman [2003] J. C. Hayes and M. L. Norman, Beyond Flux-limited Diffusion: Parallel Algorithms for Multidimensional Radiation Hydrodynamics, The Astrophysical Journal Supplement Series 147, 197 (2003), publisher: American Astronomical Society.
  * Hubeny and Burrows [2007] I. Hubeny and A. Burrows, A New Algorithm for Two-Dimensional Transport for Astrophysical Simulations. I. General Formulation and Tests for the One-Dimensional Spherical Case, The Astrophysical Journal 659, 1458 (2007), publisher: American Astronomical Society.
  * Dasgupta  _et al._ [2008] B. Dasgupta, A. Dighe, A. Mirizzi, and G. Raffelt, Collective neutrino oscillations in nonspherical geometry, Phys.Rev.D 78, 033014 (2008), arXiv:0805.3300 [hep-ph] .
  * Mirizzi and Serpico [2012] A. Mirizzi and P. D. Serpico, Instability in the Dense Supernova Neutrino Gas with Flavor-Dependent Angular Distributions, Physical Review Letters 108, 231102 (2012), publisher: American Physical Society.
  * Tamborra and Shalgar [2020] I. Tamborra and S. Shalgar, New Developments in Flavor Evolution of a Dense Neutrino Gas, arXiv e-prints , arXiv:2011.01948 (2020), arXiv:2011.01948 [astro-ph.HE] .
  * Duan  _et al._ [2007] H. Duan, G. M. Fuller, J. Carlson, and Y.-Z. Qian, Neutrino Mass Hierarchy and Stepwise Spectral Swapping of Supernova Neutrino Flavors, Physical Review Letters 99, 241802 (2007), publisher: American Physical Society.
  * Dasgupta  _et al._ [2009a] B. Dasgupta, A. Dighe, G. G. Raffelt, and A. Y. Smirnov, Multiple Spectral Splits of Supernova Neutrinos, Physical Review Letters 103, 051105 (2009a), publisher: American Physical Society.
  * Keil  _et al._ [2003] M. T. Keil, G. G. Raffelt, and H. Janka, Monte Carlo Study of Supernova Neutrino Spectra Formation, The Astrophysical Journal 590, 971 (2003).
  * Fischer  _et al._ [2010] T. Fischer, S. C. Whitehouse, A. Mezzacappa, F. K. Thielemann, and M. Liebendörfer, Protoneutron star evolution and the neutrino-driven wind in general relativistic neutrino radiation hydrodynamics simulations, Astronomy & Astrophysics 517, A80 (2010), arXiv:0908.1871 [astro-ph.HE] .
  * Zyla  _et al._ [2020] P. Zyla _et al._ (Particle Data Group), Review of Particle Physics, PTEP 2020, 083C01 (2020).
  * Duan  _et al._ [2006a] H. Duan, G. M. Fuller, and Y.-Z. Qian, Collective neutrino flavor transformation in supernovae, Physical Review D 74, 123004 (2006a), publisher: American Physical Society.
  * Dasgupta  _et al._ [2009b] B. Dasgupta, A. Dighe, G. G. Raffelt, and A. Y. Smirnov, Multiple spectral splits of supernova neutrinos, Phys. Rev. Lett. 103, 051105 (2009b).
  * Hunter [2007] J. D. Hunter, Matplotlib: A 2d graphics environment, Computing in Science Engineering 9, 90 (2007).
  * van der Walt  _et al._ [2011] S. van der Walt, S. C. Colbert, and G. Varoquaux, The numpy array: A structure for efficient numerical computation, Computing in Science Engineering 13, 22 (2011).
  * Jones  _et al._ [2001] E. Jones, T. Oliphant, P. Peterson, _et al._ , SciPy: Open source scientific tools for Python (2001), [Online; accessed ¡today¿].
  * Duan  _et al._ [2006b] H. Duan, G. M. Fuller, J. Carlson, and Y.-Z. Qian, Coherent Development of Neutrino Flavor in the Supernova Environment, Physical Review Letters 97, 241101 (2006b), publisher: American Physical Society.
  * Duan  _et al._ [2006c] H. Duan, G. M. Fuller, J. Carlson, and Y.-Z. Qian, Simulation of coherent non-linear neutrino flavor transformation in the supernova environment. i: Correlated neutrino trajectories, Phys. Rev. 74, 105014 (2006c), arXiv:astro-ph/0606616 .
  * Keister [2015] B. D. Keister, Numerical and physical stability of supernova neutrino flavor evolution, Physica Scripta 90, 088008 (2015), arXiv:1408.4729 [astro-ph.HE] .
  * Richers  _et al._ [2019] S. A. Richers, G. C. McLaughlin, J. P. Kneller, and A. Vlasenko, Neutrino quantum kinetics in compact objects, Phys.Rev.D 99, 123014 (2019), arXiv:1903.00022 [astro-ph.HE] .