# Two-Moment Neutrino Flavor Transformation with applications to the Fast
Flavor Instability in Neutron Star Mergers

Evan Grohs  Department of Physics, North Carolina State University, Raleigh,
NC 27695, USA  Sherwood Richers  Department of Physics and Astronomy,
University of Tennessee Knoxville, Knoxville, TN 37996, USA  Sean M. Couch
Department of Physics and Astronomy, Michigan State University, East Lansing,
Michigan 48824, USA  Department of Computational Mathematics, Science, and
Engineering, Michigan State University, East Lansing, Michigan 48824, USA
Facility for Rare Isotope Beams, Michigan State University, East Lansing,
Michigan 48824, USA  Francois Foucart  Department of Physics & Astronomy,
University of New Hampshire, 9 Library Way, Durham NH 03824, USA  Julien
Froustey  Department of Physics, North Carolina State University, Raleigh, NC
27695, USA  Department of Physics, University of California Berkeley,
Berkeley, CA 94720, USA  Jim Kneller  Department of Physics, North Carolina
State University, Raleigh, NC 27695, USA  Gail McLaughlin  Department of
Physics, North Carolina State University, Raleigh, NC 27695, USA

###### Abstract

Multi-messenger astrophysics has produced a wealth of data with much more to
come in the future. This enormous data set will reveal new insights into the
physics of core collapse supernovae, neutron star mergers, and many other
objects where it is actually possible, if not probable, that new physics is in
operation. To tease out different possibilities, we will need to analyze
signals from photons, neutrinos, gravitational waves, and chemical elements.
This task is made all the more difficult when it is necessary to evolve the
neutrino component of the radiation field and associated quantum-mechanical
property of flavor in order to model the astrophysical system of interest — a
numerical challenge that has not been addressed to this day. In this work, we
take a step in this direction by adopting the technique of angular-integrated
moments with a truncated tower of dynamical equations and a closure,
convolving a flavor-transformation with spatial transport to evolve the
neutrino radiation quantum field. We show that moments capture the dynamical
features of fast flavor instabilities and provide comparable results to a more
precise particle-in-cell method. We propose areas for improvement in the
future.

Neutrino flavor transformation; angular moments; neutron star mergers

N3AS-23-029

##  1 Introduction

Both neutron star mergers and core-collapse supernovae are true multi-
messenger events, as they produce neutrinos, photons, gravitational waves, and
chemical elements. In the coming decade, there will be a wealth of data from
all of these messengers, see e.g. Kalogera et al. (2021); Baxter et al.
(2022); Holmbeck et al. (2020); Cowperthwaite et al. (2017); Lien & Fields
(2009); LSST Science Collaboration et al. (2009); Bellm et al. (2018);
Tartaglia et al. (2018). In order to produce the most realistic theoretical
predictions to compare with future data, much theoretical development is still
needed. One significant area that requires attention is the neutrino physics
of hot and dense systems (for a recent review see Volpe 2023).

Stellar explosions that reach extreme temperatures and densities, such as
core-collapse supernovae and neutron star mergers, produce enough neutrinos
that they account for a substantial portion of the energy budget (for recent
estimates see Bollig et al. 2021; Burrows et al. 2019; Fujibayashi et al.
2023; Hayashi et al. 2022; Foucart et al. 2023; Cusinato et al. 2022). The
majority of these neutrinos are in the energy range of tens of MeV. In
neutrino rich regions, the ratio of neutrons to protons is influenced by
electron neutrino and electron antineutrino capture reactions. This neutron-
to-proton ratio is a key factor influencing element synthesis (e.g.,
McLaughlin et al. 1996; Freiburghaus et al. 1999; Surman et al. 2006; Lippuner
& Roberts 2015; Curtis et al. 2019; Miller et al. 2020; Reichert et al. 2021;
Curtis et al. 2023).

Exploratory work has demonstrated the importance of accurately understanding
the impact of changes in neutrino flavor as a function of both time and
position in the exploding object. For example, the distribution of neutrinos
among flavors influences the outcome of the supernova explosions in one-
dimensional high-mass supernovae simulations (Stapleford et al., 2020; Ehring
et al., 2023a, b), the outcome of element synthesis obtained in neutrino-
cooled accretion disks (e.g., Malkus et al. 2012; Just et al. 2022), in
supernovae (e.g., Duan et al. 2011; Mukhopadhyay 2022) and in hypermassive
neutron star outflows (e.g., Fernández et al. 2022). Since the type of
elements that are produced in ejecta depend sensitively on the ratio of
neutrons to protons, these studies conclude that there is an impact on the
elements that are produced in Core-Collapse SuperNovae (CCSNe) and Neutron
Star Mergers (NSMs).

The Quantum Kinetic Equations (QKEs), where the terms representing the
interactions of the neutrinos are expanded in a series, are often taken as a
starting point for calculating the outcome of neutrino transport and
propagation. The first term in this series corresponds to evolution through a
potential while the second term corresponds to momentum-changing collisions
(e.g., Volpe et al. 2013; Vlasenko et al. 2014; Blaschke & Cirigliano 2016;
Froustey et al. 2020). Different groups use the phrase QKE in a variety of
manners, referring to the specific terms included/excluded in the series
expansion. In this work, we will use the phrase QKE to denote any equation
which modifies the neutrino density matrices in time. The starting point for
classical neutrino transport can be obtained from the neutrino QKEs under the
approximation that neutrino density matrices are always on-diagonal. In this
case, only the second term in the QKE series expansion – the collision term –
is relevant.

Modern codes aiming to perform global 3D simulations of CCSNe (Just et al.,
2015; Kuroda et al., 2016; Skinner et al., 2019; Bruenn et al., 2020) and/or
NSMs and post-merger remnants (Ruffert & Janka, 1999; Rosswog & Liebendörfer,
2003; Wanajo et al., 2014; Neilsen et al., 2014; Perego et al., 2016; Foucart
et al., 2016a; Ardevol-Pulpillo et al., 2019; Gizzi et al., 2019; Foucart et
al., 2021; Radice et al., 2022) that incorporate neutrinos generally use such
classical transport algorithms. Given the difficulty of accurately solving the
transport problem in global simulations with sufficient resolution and
detailed microphysics, these codes inevitably need to make a range of
additional approximations. On the methods side, these might include the use of
approximate transport schemes (leakage, truncated moments) or low-resolution
Monte-Carlo methods, and/or the use of energy-integrated transport or ray-by-
ray transport. On the microphysics side, this often includes the use of
approximate interaction rates, reduced number of neutrino flavors, or simply
ignoring interactions that are too costly to calculate in practice, e.g., pair
processes and inelastic scattering. While global neutrino transport algorithms
are rapidly improving, they are still having significant difficulties in
capturing all important aspects of the classical transport equation.

Meanwhile, the starting point for flavor transformation in the absence of
collisions is often studied by evolving the flavor field using only the first
term in the QKE series expansion. The evaluation of this term is done by use
of operator splitting of the Hamiltonian (so-called _mean-field_). This set-up
has been studied extensively: for example, the part of this Hamiltonian
associated with neutrino coherent forward scattering on other neutrinos, in
combination with other Hamiltonian terms, gives rise to the phenomenon of
bipolar oscillations (for a review see Duan et al. 2006) and matter neutrino
resonance transitions (Malkus et al., 2012, 2014; Frensel et al., 2017; Wu et
al., 2016; Tian et al., 2017; Vlasenko & McLaughlin, 2018). Additionally, Fast
Flavor Conversion (FFC) stems from the combination of specific angular
distributions of neutrinos, the mean-field Hamiltonian, and inclusion of
neutrino advection, e.g. Sawyer (2005); Dasgupta et al. (2017); Izaguirre et
al. (2017). The relevant angular distributions are expected to occur in both
supernovae, e.g. Abbar et al. (2019); Nagakura et al. (2021); Nagakura (2023),
and neutron star mergers (Wu & Tamborra, 2017) at positions close to the
central object. A number of works exist that evaluate classically computed
angular distributions to determine whether these distributions have a Fast
Flavor Instability (FFI), using a variety of techniques (Richers, 2022; Abbar,
2023). The hallmark of a test for whether an instability will exist is to look
for a “crossing” between a curve that represents the number density of
neutrinos as a function of angle, and the curve that represents the number
density of antineutrinos as a function of angle, e.g. Dasgupta et al. (2009);
Abbar & Duan (2018); Dasgupta (2022).

Ideally one wishes to use both of the first two terms in the quantum kinetic
equation series and some efforts have been undertaken with the inclusion of
both. When including both the first and second term in the QKE series, the
collision term most typically produces decoherence of the neutrinos, e.g.
Richers et al. (2019), and if that term is sufficiently large, the neutrinos
tend to drift into flavor states. However, under the right conditions, the
combination of the two terms can also enhance flavor transformation through
collisional instabilities (Johns, 2023; Johns & Xiong, 2022; Xiong et al.,
2023, 2022).

At present, there are questions about whether the QKEs can ever completely
capture the behavior of neutrinos in these astrophysical systems, specifically
because of the operator splitting technique that is used to write down the
Hamiltonian. There are ongoing efforts to analyze the evolution of neutrinos
due to the forward scattering part of the many-body Hamiltonian under the
assumption of continuous temporal interactions of all neutrinos with all other
neutrinos (Balantekin & Pehlivan, 2007; Pehlivan et al., 2011; Siwach et al.,
2023; Lacroix et al., 2022; Balantekin et al., 2023; Cervia et al., 2022;
Rrapaj, 2020; Patwardhan et al., 2021; Roggero et al., 2022; Martin et al.,
2023b, a).

Notwithstanding these questions, efforts have been made to compute the QKEs by
capturing the evolution of many neutrino “packets” in many different
directions (e.g., Sawyer 2005; Martin et al. 2020; Richers et al. 2021a;
Bhattacharyya & Dasgupta 2020; Zaizen & Morinaga 2021; Nagakura 2022; George
et al. 2022). These methods provide useful benchmarks but are at present too
computationally expensive to use extensively.

An alternative is to use a reformulation of the QKEs in terms of the angular
moments (Zhang & Burrows, 2013; Richers et al., 2019). This reformulation
creates a series of equations describing the time evolution of each moment and
one then evolves only a small number of these equations. One then has to
choose what to do with the moments which appear in the evolution equations but
which are not explicitly evolved. One approach is simply to ignore the
evolution of the moments above some order, however, what is found in practice
is that one needs to retain a large number of the moment evolution equations
making this approach computationally inefficient (Dasgupta et al., 2018; Johns
et al., 2020, 2020). An alternative solution to the truncation problem is to
use a “closure” which links the unevolved moments to the lower order, evolved
moments via some relationship. An example of calculations using this closure
method can be found in Myers et al. (2022), and using this approach, moment
methods have been able to reproduce fast flavor transformations in neutron
star merger-like conditions (Grohs et al., 2022).

In this manuscript, we consider in detail the efficacy of a two moment
implementation of the QKEs neglecting the collision term and using a closure.
We illustrate our method using an example quantum closure that is a relatively
straightforward extension of the classical maximum entropy closure. In Sec. 2
we detail the QKE formalism and apply it to moments. We elucidate the angular
distributions corresponding to the closure and how they imply lepton number
crossings in Sec. 3. Section 4 gives an exposition of our implementation of
neutrino flavor transformation in the framework of FLASH (Fryxell et al.,
2000; Dubey et al., 2009) along with the initial and boundary conditions. In
Sec. 5 we compare the results of our moment treatment to a more exact method
for several well-studied test problems before turning to the presentation of
results for three kinds of neutron star merger-like conditions. We conclude
and discuss the need for further exploration of moment QKE methods in Sec. 6.
With regard to units: we use two conventions. When writing the neutrino flavor
transformation equations, we use natural units where $\hbar=c=1$. When giving
results of numerical calculations, we use cgs units.

##  2 The Moment Evolution Equations

###  2.1 General Formalism

We begin from the general QKEs describing the neutrino and anti-neutrino
evolution adopting the nomenclature of Sigl & Raffelt (1993); Vlasenko et al.
(2014); Froustey et al. (2020) and in particular Blaschke & Cirigliano (2016).
The evolved variable in the QKEs is a generalized density matrix for
neutrinos, $\varrho=\varrho(t,\mathbf{x},\mathbf{p})$, and corresponding
generalized density matrix $\overline{\varrho}$ for anti-neutrinos, which are
functions of time $t$, spatial location $\mathbf{x}$, and momentum
$\mathbf{p}$. In the treatment of Blaschke & Cirigliano (2016), the
generalized density matrices are one-body reduced density matrices (Volpe et
al., 2013; Froustey et al., 2020), and hereafter we will call $\varrho$ and
$\overline{\varrho}$ simply “density matrices” for the sake of brevity. If we
are describing neutrinos and anti-neutrinos with 2 chiral states, $\varrho$
and $\overline{\varrho}$ are $2\,n_{f}\times 2\,n_{f}$ Hermitian matrices for
$n_{f}$ flavors. However in this work, we only consider left-chiral neutrinos
and right-chiral anti-neutrinos, and ignore any kind of spin coherence
(Cirigliano et al., 2015; Tian et al., 2017). As a result, the size of the
density matrices is reduced to $n_{f}\times n_{f}$ for each of the neutrinos
and anti-neutrinos. In the case of three flavors of neutrinos, namely
$e,\mu,\tau$, we can write the neutrino density matrix as the following

$\varrho=\begin{pmatrix}\varrho_{ee}&\varrho_{e\mu}&\varrho_{e\tau}\\\
\varrho_{\mu e}&\varrho_{\mu\mu}&\varrho_{\mu\tau}\\\ \varrho_{\tau
e}&\varrho_{\tau\mu}&\varrho_{\tau\tau},\end{pmatrix}$ (1)

with a similar expression for anti-neutrinos. In Eq. (1), the diagonal terms
indicate the occupation numbers for a given flavor. The off-diagonal terms of
Eq. (1) encode the quantum coherence between two flavors. The expressions for
the number density, energy density, and the number flux vector are obtained by
taking appropriate phase-space integrals of $\varrho$:

$\displaystyle\mathcal{N}(t,\mathbf{x})$
$\displaystyle=\frac{1}{(2\pi)^{3}}\int
d^{3}p\,\varrho(t,\mathbf{x},\mathbf{p})\,,$ (2)
$\displaystyle\mathcal{E}(t,\mathbf{x})$
$\displaystyle=\frac{1}{(2\pi)^{3}}\int
d^{3}p\,p\,\varrho(t,\mathbf{x},\mathbf{p})\,,$ (3)
$\displaystyle\mathcal{F}^{i}(t,\mathbf{x})$
$\displaystyle=\frac{1}{(2\pi)^{3}}\int
d^{3}p\,\frac{p^{i}}{p}\,\varrho(t,\mathbf{x},\mathbf{p})\,,$ (4)

where the superscript $i$ indicates a component of a 3-vector. Note that we
approximate neutrinos as ultra-relativistic by setting the neutrino energy
equal to the 3-momentum magnitude $p$ and that $\mathcal{N}$, $\mathcal{E}$,
and the components of $\mathbf{\mathcal{F}}$ are all $n_{f}\times n_{f}$
matrices. The expressions for anti-neutrinos are analogous.

The QKEs for $\varrho$ and $\overline{\varrho}$ are

$\displaystyle\frac{\partial\varrho}{\partial
t}+\dot{\mathbf{x}}\cdot\frac{\partial\varrho}{\partial\mathbf{x}}+\dot{\mathbf{p}}\cdot\frac{\partial\varrho}{\partial\mathbf{p}}$
$\displaystyle=-\imath\,[H,\varrho]+C,$ (5)
$\displaystyle\frac{\partial\overline{\varrho}}{\partial
t}+\dot{\mathbf{x}}\cdot\frac{\partial\overline{\varrho}}{\partial\mathbf{x}}+\dot{\mathbf{p}}\cdot\frac{\partial\overline{\varrho}}{\partial\mathbf{p}}$
$\displaystyle=-\imath\,[\overline{H},\overline{\varrho}]+\overline{C},$ (6)

where the single dot over a variable indicates differentiation with respect to
time. In Eqs. (5) and (6), $C$ and $\overline{C}$ are collision terms which
can change neutrino flavor, number, or momenta. We shall ignore them
throughout this work given the large separation of scales between the fast-
flavor instability growth rate and the collision rates simulated here. In
addition, we will consider systems where the particle 3-momenta do not change
with time, implying we may exclude the force term on the lhs of Eqs. (5) and
(6). To study flavor transformation, we employ Hamiltonian-like operators in
Eqs. (5) and (6) consistent with mean-field treatments. When working to first
order in power counting of the QKEs (Vlasenko et al., 2014), the Hamiltonian
operators are a sum of three potentials. Specifically,

$\displaystyle H$ $\displaystyle=H_{V}+H_{M}+H_{\nu},$ (7)
$\displaystyle\overline{H}$ $\displaystyle=H_{V}-H_{M}-H_{\nu}^{\ast},$ (8)

denoting the vacuum $(H_{V})$, matter $(H_{M})$, and self-interaction
$(H_{\nu})$ terms, and where ∗ denotes complex conjugation. The vacuum term
arises from non-zero neutrino rest masses, and we write it as

$H_{V}=\frac{1}{2p}\,U\,M^{2}\,U^{\dagger}\,,$ (9)

where $U$ is the PMNS matrix and
$M^{2}=\mathrm{diag}(m_{1}^{2},m_{2}^{2},m_{3}^{2})$ is the diagonal matrix of
squared neutrino masses. The matter term is linear and familiar in the context
of oscillations with solar neutrinos. Electrons and positrons interact weakly
with neutrinos in a flavor-dependent manner, which we denote by the following
expression

$H_{M}=\sqrt{2}G_{F}\,n_{e}\,I_{e},$ (10)

where $G_{F}\simeq 1.166\times 10^{-11}\,{\rm MeV}^{-2}$ is the Fermi coupling
constant, $n_{e}$ is the difference between the number density of electrons
and positrons, and $I_{e}$ is the electron flavor projection operator, i.e.
$I_{e}={\rm diag}(1,0,0)$ for three flavors. Finally, the self-interaction
potential is a consequence of neutrinos interacting with the background of
other neutrinos

$H_{\nu}=\frac{\sqrt{2}G_{F}}{(2\pi)^{3}}\int
d^{3}q(1-\cos\vartheta)[\varrho(t,\mathbf{x},\mathbf{q})-\overline{\varrho}^{\ast}(t,\mathbf{x},\mathbf{q})],$
(11)

where $\vartheta$ is the angle between the free variable $\mathbf{p}$ and the
integration variable $\mathbf{q}$.

###  2.2 Moment Quantum Kinetic Equations

In general, the density matrices are seven-dimensional since they depend upon
time, space, and momentum. Solving the QKEs for the density matrices with
sufficient temporal, spatial, and momentum resolution to ensure numerical
convergence will be very computationally expensive. An alternative approach is
to recast the QKEs as an infinite set of transport equations for the moments
of the density matrices, and then truncate the number of moments that one
solves by adopting a closure. Since moments are only five-dimensional
quantities, solving their transport equations with sufficient fidelity to
ensure convergence is a more feasible, though still difficult, computational
challenge. In this paper we adopt a two moment scheme in which we evolve only
the “zeroth” and “first” angular-integrated moments. We define the zeroth,
first and second moment of $\varrho$ to be

$\displaystyle E(t,\mathbf{x},p)$
$\displaystyle=\frac{1}{4\pi}\frac{p^{3}}{(2\pi)^{3}}\int
d\Omega_{p}\,\varrho(t,\mathbf{x},\mathbf{p}),$ (12) $\displaystyle
F^{i}(t,\mathbf{x},p)$
$\displaystyle=\frac{1}{4\pi}\frac{p^{3}}{(2\pi)^{3}}\int
d\Omega_{p}\frac{p^{i}}{p}\,\varrho(t,\mathbf{x},\mathbf{p}),$ (13)
$\displaystyle P^{ij}(t,\mathbf{x},p)$
$\displaystyle=\frac{1}{4\pi}\frac{p^{3}}{(2\pi)^{3}}\int
d\Omega_{p}\frac{p^{i}p^{j}}{p^{2}}\,\varrho(t,\mathbf{x},\mathbf{p}),$ (14)

where $i,j\in\\{x,y,z\\}$ are spatial indices. Note that the factor of
$1/4\pi$ which appears in the definitions of these moments is a units choice
made to divide out the $4\pi$ from the momentum-space solid angle integral.
Analogous expressions exist for the anti-neutrinos. The integrals in these
definitions are only over the momentum-space solid angle $\Omega_{p}$ i.e. the
propagation directions of the neutrinos at a given spacetime location and
comoving-frame neutrino energy, and not the entire phase-space as in Eqs. (2)
– (4). Furthermore, the $p^{3}$ in the prefactor of Eqs. (12) – (14) indicates
that these are the differential energy density and differential energy flux.
These are the quantities we have chosen to time-evolve in the FLASH code since
this is the convention used in many instances of classical neutrino moment
transport. However we will oftentimes show results using instead the
differential number density which is related to $E(t,\mathbf{x},p)$ via

$N(t,\mathbf{x},p)=\frac{1}{p}\,E(t,\mathbf{x},p).$ (15)

From the definition of the moments we see we can recover the expressions in
Eqs. (2) – (4) by integrating the moments over the neutrino energy $p$ i.e.

$\displaystyle\mathcal{N}(t,\mathbf{x})$ $\displaystyle=4\,\pi\,\int
dp\,N(t,\mathbf{x},p)\,,$ (16) $\displaystyle\mathcal{E}(t,\mathbf{x})$
$\displaystyle=4\,\pi\,\int dp\,E(t,\mathbf{x},p)\,,$ (17)
$\displaystyle\mathcal{F}^{i}(t,\mathbf{x})$ $\displaystyle=4\,\pi\,\int
dp\,\frac{F^{i}(t,\mathbf{x},p)}{p}\,.$ (18)

Note that $F$ is the specific energy flux, but $\mathcal{F}$ is the energy-
integrated number flux. For future reference, at this point we introduce the
flux factor vector (actually a vector of matrices) which we define to be

$\mathbf{f}_{ab}=\frac{\mathbf{F}_{ab}}{E_{ab}},$ (19)

with the norm for a component of the flavor matrix defined as

$f_{ab}\equiv|\mathbf{f}_{ab}|=\sqrt{\sum\nolimits_{i}(f^{i}_{ab})^{2}}.$ (20)

From comparing the definitions of the moments and Eq. (11) for self-
interactions, we observe that the self-interaction term can be written as

$H_{\nu}=H_{E}-\frac{1}{p}\,\mathbf{p}\cdot\mathbf{H}_{F},$ (21)

where the moment-self-interaction terms are

$\displaystyle H_{E}$
$\displaystyle=\sqrt{2}G_{F}\left(\mathcal{N}-\overline{\mathcal{N}}^{*}\right),$
(22) $\displaystyle H_{F}^{j}$
$\displaystyle=\sqrt{2}G_{F}\left(\mathcal{F}^{j}-\overline{\mathcal{F}}^{*j}\right).$
(23)

Although the physical systems which we model do occur in environments where
general relativity has a pronounced effect over large distances, we will do
all of our calculations in a local Minkowski reference frame where we may
specify the spacetime metric as $g^{\mu\nu}={\rm diag}(-1,1,1,1)$. Therefore,
the 3-vector contraction in Eq. (21) (and all other subsequent 3-vector
contractions in this work) is equivalent to the 3D dot product of Euclidean
space.

With everything defined we can now write out the evolution equations for $E$
and $\mathbf{F}$ by performing moment integrations of Eqs. (5) and (6) and
scaling by appropriate factors of $p^{3}/(2\pi)^{3}/(4\pi)$. In Cartesian
coordinates the transport equations for the neutrino moments are

$\displaystyle\frac{\partial E}{\partial t}+\frac{\partial F^{j}}{\partial
x^{j}}$
$\displaystyle=-\imath\,[H_{V}+H_{M}+H_{E},E]+\imath\,[H_{F}^{j},F_{j}],$ (24)
$\displaystyle\frac{\partial F^{j}}{\partial t}+\frac{\partial
P^{jk}}{\partial x^{k}}$
$\displaystyle=-\imath\,[H_{V}+H_{M}+H_{E},F^{j}]+\imath\,[H_{F}^{k},P^{j}_{k}],$
(25)

where we have ignored the collision and force terms. Note that the second
(pressure) moment is not an evolved quantity in a M1 transport scheme and is
calculated algebraically as a function of the two time-evolved moments $E$ and
$\mathbf{F}$ using a closure relation. The closure relation we use will be
discussed in Sec. 3.1. Finally for completeness, we give the moment evolution
equations for the anti-neutrinos

$\displaystyle\frac{\partial\overline{E}}{\partial
t}+\frac{\partial\overline{F}^{j}}{\partial x^{j}}$
$\displaystyle=-\imath\,[H_{V}-H_{M}-H_{E}^{\ast},\overline{E}]+\imath\,[H_{F}^{j\ast},\overline{F}_{j}],$
(26) $\displaystyle\frac{\partial\overline{F}^{j}}{\partial
t}+\frac{\partial\overline{P}^{jk}}{\partial x^{k}}$
$\displaystyle=-\imath\,[H_{V}-H_{M}-H_{E}^{\ast},\overline{F}^{j}]+\imath\,[H_{F}^{k\ast},\overline{P}^{j}_{k}].$
(27)

Equations (24), (25), (26), and (27) give the equations of motion for the
neutrino field under study. They comprise a coupled set of matrix equations
with spacetime indices $\\{j,k\\}$ ranging over the 3D space indices
$\\{1,2,3\\}$. As an interesting aside, we note that it is possible to write
the equations of motion in a 4D spacetime framework. To begin, construct the
following neutrino arrays

$\displaystyle J^{\alpha}_{(\nu)}$
$\displaystyle=\begin{pmatrix}\mathcal{N}-\overline{\mathcal{N}}^{*}\\\
\mathcal{F}-\overline{\mathbf{\mathcal{F}}}^{*}\end{pmatrix},$ (28)
$\displaystyle T^{\alpha\beta}$ $\displaystyle=\begin{pmatrix}E&\mathbf{F}\\\
\mathbf{F}^{T}&\mathbf{P}\end{pmatrix},$ (29) $\displaystyle H^{\alpha}$
$\displaystyle=-u^{\alpha}H_{V}-\sqrt{2}G_{F}(J^{\alpha}_{(e)}+J^{\alpha}_{(\nu)}),$
(30) $\displaystyle\overline{H}^{\alpha}$
$\displaystyle=-u^{\alpha}H_{V}+\sqrt{2}G_{F}(J^{\alpha}_{(e)}+J^{\alpha}_{(\nu)})^{*},$
(31)

where $J^{\alpha}_{(e)}=u^{\alpha}(n_{e}-\bar{n}_{e})I_{e}$. In the above Eq.
(29), $\mathbf{F}^{T}$ denotes the transpose of the row vector $\mathbf{F}$
into a column vector. In addition, we define the four-velocity of the
reference frame $u^{\mu}=(1,0,0,0)$ in Eqs. (30) and (31). With these
definitions, we are able to cast the QKEs for neutrinos and anti-neutrinos as

$\displaystyle\nabla_{\beta}T^{\alpha\beta}$
$\displaystyle=-i\left[H_{\beta},T^{\alpha\beta}\right],$ (32)
$\displaystyle\nabla_{\beta}\overline{T}^{\alpha\beta}$
$\displaystyle=-i\left[\overline{H}_{\beta},\overline{T}^{\alpha\beta}\right].$
(33)

In Eqs. (32) and (33), we adopt the convention where repeated indices are
contracted with respect to the metric, i.e.,
$A^{\alpha}A_{\alpha}=g_{\alpha\beta}A^{\alpha}A^{\beta}$ with a $(-,+,+,+)$
convention.

##  3 Moment Closure Relation and Lepton Number Crossing

###  3.1 The Maximum Entropy Closure

The evolution equations for the fluxes $\mathbf{F}$ and
$\overline{\mathbf{F}}$ involve the spatial gradients of the pressure tensors
$P$ and $\overline{P}$; the evolution of the pressure tensors involve the
spatial gradients of the next moment. This pattern continues in perpetuity and
results in an infinite tower of equations. This is an unavoidable property of
moment decomposition. Nevertheless, in some situations the infinite set of
equations can be solved: for example, when the radiation field is strongly-
interacting, an equation of state will relate the pressure to the energy
density under the assumption of Local Thermodynamic Equilibrium (LTE) thereby
closing the set of equations for the first two moments. But in general – and
neutrinos in CCSNe and NSMs are both such cases – no such equation of state
exists that naturally closes the set of evolution equations. The simplest
approach is to propose a local, analytic relation to close the tower of
equations suited for the individual problem under study that matches analytic
results in the trapped and free-streaming limits. This relation is called the
closure relation (or “closure” for brevity). We will adopt this same approach
when proposing a closure for the quantum moments of the neutrinos that must be
able to account for both neutrino advection and the flavor transformation. We
begin our explanation of the closure we adopt by ignoring the flavor structure
of the moments for the time being, and consider the Maximum Entropy Closure
(MEC) often used in “classical” moment transport.

By definition of the MEC, the neutrinos of a particular species assume an
angular distribution in momentum-space such that the angular entropy is
extremized (Minerbo, 1978; Cernohorsky & Bludman, 1994). In other words, the
neutrinos are distributed in the momentum-space angles such that an entropy-
like function is maximized under the constraints of a net number density and
flux. These constraints relate directly to the dynamical variables of interest
in Eqs. (12) and (13). As with any reasonable closure, the MEC exactly
represents the radiation field in the limits far from a source (where all
radiation is moving in one direction) and when the radiation is in
equilibrium. We will utilize the MEC for our flavor-mixing neutrino and anti-
neutrino distributions, which we summarize below for completeness.

Under the constraints of number density and flux, the neutrino distribution of
a particular species $a$ per unit solid angle of momentum-space, $\psi_{aa}$,
is (Minerbo, 1978; Cernohorsky & Bludman, 1994)

$\psi_{aa}=\frac{E_{aa}}{4\pi}\,\frac{Z_{aa}}{\mathrm{sinh}(Z_{aa})}e^{Z_{aa}\mu}\,\,,$
(34)

where $E_{aa}$ is the energy density moment for species $a$, and
$\mu=\hat{F}_{\mathrm{tet}}\cdot\mathbf{\Omega}$ gives the angular dependence.
$\mathbf{\Omega}$ is the direction unit vector and $\widehat{F}_{\rm
tet}=\mathbf{F}_{aa}/F_{aa}$. The parameter $Z_{aa}$ follows from the
constraint on the magnitude of the flux factor vector $f_{aa}$

$\displaystyle f_{aa}$ $\displaystyle=\frac{1}{E_{aa}}\int
d\Omega\,\mu\,\psi_{aa}$ (35) $\displaystyle=\coth
Z_{aa}-\frac{1}{Z_{aa}}\,\,.$

Eq. (35) must, in general, be inverted numerically to obtain the value of
$Z_{aa}$ corresponding to a given $f_{aa}$. Once $Z_{aa}$ is obtained, we can
construct the angular distribution of the neutrinos in Eq. (34). If we adopt
initial neutrino distributions from a core collapse supernova or neutron star
merger simulation that uses an MEC, then this is consistent with the
assumptions in the original simulations. While the full angular information is
assumed in our multi-direction calculations, our two-moment scheme simply uses
the MEC to determine the pressure moment.

Borrowing the terminology from classical radiation hydrodynamics, we
interpolate the pressure moment between the optically thin and thick limits as

$P^{ij}=\frac{3\chi-1}{2}P^{ij}_{\mathrm{thin}}+\frac{3(1-\chi)}{2}P^{ij}_{\mathrm{thick}}\,\,,$
(36)

where the thin and thick limits are

$\displaystyle P_{\mathrm{thin}}^{ij}$
$\displaystyle=E\frac{F^{i}F^{j}}{F^{2}},$ (37) $\displaystyle
P_{\mathrm{thick}}^{ij}$ $\displaystyle=\frac{E}{3}\delta^{ij},$ (38)

and we have suppressed the flavor indices for ease in notation. (Minerbo,
1978) demonstrate that these assumptions lead to a simple functional form of
the Eddington factor $\chi$:

$\chi=\frac{1}{3}+\frac{2}{15}f^{2}(3-f+3f^{2}),$ (39)

such that Eq. (36) becomes consistent with the second angular moment of Eq.
(34). In Sec. 4.1 we explore extending this concept to matrix-valued moments
necessary for quantum neutrino transport.

###  3.2 Lepton Number Crossing with the Maximum Entropy Closure

At this point, we give a brief interlude to discuss lepton number crossings in
the context of the MEC. Assuming the momentum-space angular distributions of
two neutrino species follow an MEC, one can analytically determine whether two
distributions cross (Johns & Nagakura, 2021; Richers, 2022). Such crossings
(most straightforwardly between electron neutrino and anti-neutrino
distributions) herald in flavor instability Morinaga (2022). Although the
conditions for neutrino flavor instability are more general and involve the
other flavor-lepton-numbers, we shall consider initial conditions where the
$x$-flavor Lepton Number (XLN) is zero and therefore the ELN crossing is the
sole source of the instability.

The intersection of the angular distributions is the boundary of a 2D surface
in the 3D momentum-space. Solving for the boundary is, in general a difficult
problem. However, for the purposes of FFC, simply identifying whether or not
the intersection exists suffices to determine whether the system is unstable
or not. Therefore, we can look at the 2D cross-sectional slice of the 3D
distributions in the plane of both flux vectors to determine whether there is
a ELN crossing or not.

Let $E_{ee}$ and $Z_{ee}$ define the maximum entropy distributions for
electron neutrinos, and similarly $\overline{E}_{ee}$ and $\overline{Z}_{ee}$
for electron anti-neutrinos [see Eqs. (34) and (35)]. Furthermore, we assume
that the flux factors are separated by an angle $\theta$, i.e.,

$\cos\theta=\frac{\mathbf{f}_{ee}\cdot\overline{\mathbf{f}}_{ee}}{f_{ee}\overline{f}_{ee}}$
(40)

The distributions cross if (Richers, 2022)

$\frac{\eta^{2}}{\alpha^{2}+\gamma^{2}}\leq 1\,\,,$ (41)

where

$\displaystyle\eta$
$\displaystyle=\ln\left[\frac{E_{ee}Z_{ee}\sinh(\overline{Z}_{ee})}{\overline{E}_{ee}\overline{Z}_{ee}\sinh(Z_{ee})}\right],$
(42) $\displaystyle\alpha$ $\displaystyle=\bar{Z}_{ee}\sin\theta,$ (43)
$\displaystyle\gamma$ $\displaystyle=\bar{Z}_{ee}\cos\theta-Z_{ee}.$ (44)

We use the criterion in Eq. (41) to indicate the presence of FFI when choosing
the locations from NSM simulations to consider in Sec. 5.2.

##  4 Methods

We have written four QKEs (one energy density and three components for the
flux density) in Eqs. (24) and (25). Along with the equations for the anti-
neutrinos, this set of coupled matrix equations comprises 32 evolution
variables per energy bin per spatial cell. Our goal will be to integrate these
equations under the conditions of FFI to see if this method can capture the
behavior of FFC. Before presenting results of test and NSM simulations, we
give some more of the pertinent details on the numerical implementation of the
moment method into FLASH. In addition, we give a brief exposition on the
Particle-In-Cell (PIC) method in EMU and how it was tailored to compare with
FLASH.

###  4.1 FLASH

To study neutrino flavor transformation with moments, we use the FLASH
radiation hydrodynamics code (Fryxell et al., 2000; Dubey et al., 2009),
further modified by O’Connor & Couch (2018) which includes an M1 moment scheme
for classical neutrino transport. It evolves the energy and flux density
moments for three species: $\nu_{e}$, $\overline{\nu}_{e}$, and
$\nu_{\mathrm{other}}$ for all other neutrinos. We modify the classical code
by by distinguishing ${\overline{\nu}}_{\mathrm{other}}$ from
$\nu_{\mathrm{other}}$, and adding the flavor off-diagonal components of the
moments. This yields a total of eight effective species that follow from the
generalized density matrices of Sec. 2.1, and specifically a 2-flavor version
of Eq. (1). For example, these eight species for the energy density moments
are $E_{ee}$, $E_{\mu\mu}$, $\mathrm{Re}(E_{e\mu})$, and
$\mathrm{Im}(E_{e\mu})$ and the four charge-conjugate counterparts for the
anti-neutrinos. In reality, neutrinos would oscillate between $e$ and the
other two flavors, namely, $\mu$ and $\tau$. Our implementation of flavor
mixing is 2-flavor for simplicity, which artificially assumes that half of the
heavy-lepton neutrinos and antineutrinos do not participate in flavor
conversion, but the number of flavors does not alter our qualitative
conclusions and it will be possible to implement an 18-species framework for
3-flavors in the future.

We decompose the domain into $16^{3}$ blocks to parallelize the computation
over processors. We use 6 ghost cells in each dimension so communication
occurs only at the end of each full time step, given a stencil size of 2 in
each direction and 3 substeps within each full step111The stencil size is the
number of grid cells referenced when evaluating a numerical derivative. The
three-step integrator was originally designed to ensure consistency in the
hydrodynamic evolution in FLASH using a general tabulated equation of state.
Our calculations do not evolve the hydrodynamics and add an unnecessary
computational cost, but we leave the structure in place to ensure future
consistency with the full FLASH framework used for ab-initio compact object
simulations.

We extend the 3-species transport subroutines from Appendix B of O’Connor &
Couch (2018) to the 8 species needed for flavor mixing. We use the same
Harten-Lax-van Leer-Einfeldt (HLLE) Riemann solver (Harten et al., 1983) to
compute fluxes between grid cells for all 8 species. We use a first-order
method to reconstruct the interface flux and pressure values, instead of the
second order TVD reconstruction employed in O’Connor & Couch (2018). Our
advection time step is set to $0.4$ times the grid cell light crossing time.

The advection and mixing evolution is done using an operator-split method,
where the mixing derivatives are given by the rhs of Eqs. (24), (25), (26),
and (27). To calculate the commutators of $2\times 2$ matrices, we decompose
the density matrix into components as detailed above and use the commutation
relations of the Pauli matrices. Mixing is only treated locally, with the
Hamiltonian-like terms specified at a given $\mathbf{x}$, or equivalently a
given cell. Unlike the solver for the advection, we use an adaptive, explicit
5th order Runge Kutta Cash-Karp (RKCK) method (Press et al., 1993) in the
mixing subroutine closely following the implementation in Grohs et al. (2016).
The time step is determined by requiring that the difference between the
embedded 4th and 5th order solutions is smaller than one part in $10^{6}$ for
each time step and violations in unitarity (i.e., particle number
conservation) are smaller than one part in $10^{3}$.

Finally, we discuss the MEC as implemented for mixing. The Eddington factor in
Eq. (39) is a derived result from the assumptions of a classical distribution
maximizing angular entropy. Although our evolved quantities are expressed in a
particular flavor basis, the physical evolution should not be basis dependent.
Naively evaluating flux factors and Eddington factors using Eq. (39) would
break basis independence. We could diagonalize the energy density moment such
that the off-diagonal components of $E$ are zero. We would also need to apply
the same unitary transformation to each vector component of $\mathbf{F}$, but
there is no guarantee that $E$ and $F^{i}$ are simultaneously diagonalizable.
In addition, the flux factors of the flavor off-diagonal quantities are in
general complex and can be arbitrarily large or small irrespective of whether
the radiation is in the trapped or free-streaming regime, making naive flux
factors for flavor off-diagonal components a poor choice for interpolating
between these regimes.

To ameliorate these issues, we can make the assumed pressure tensor
independent of the flavor basis if we calculate a single $\chi$ for neutrinos
and a single $\overline{\chi}$ for antineutrinos using flavor-traced flux
factors. Specifically, those flavor-traced flux factors are defined as

$\displaystyle f^{(FT)}$ $\displaystyle=\frac{|{\rm Tr}[\mathbf{F}]|}{{\rm
Tr}[E]}$ (45)
$\displaystyle=\frac{|(F^{i}_{ee}+F^{i}_{xx})\hat{x}_{i}|}{E_{ee}+E_{xx}},$
(46)

and a similar expression for the anti-neutrinos and $\overline{f}^{(FT)}$,
where $\hat{x}_{i}$ are the Cartesian unit vectors. These flavor traced flux
factors are substituted into Eq. (39) to obtain $\chi$ and $\bar{\chi}$, which
are in turn used for all flavor components in Eq. (36). This also prevents the
flavor off-diagonal components from appearing to be in the optically thick
regime when the flavor off-diagonal components are in the free-streaming
regime. Note however, the principal direction of the pressure tensor is
computed as in Eq. (36) separately for _each_ flavor component, i.e.,
$\varrho_{ee}$, $\varrho_{xx}$, ${\rm Re}[\varrho_{ex}]$, and ${\rm
Im}[\varrho_{ex}]$. In other words, Eqs. (36), (37), and (38) all have flavor
indices on each quantity, except for $\chi$.

This scheme has the disadvantage that in the limit of no flavor mixing, it
does not reduce to the original two moment transport scheme, since different
flavors are no longer allowed to have independent flux factors. However, for
many of the cases we study in Sec. 5, we are in the optically thick limit for
both $e$ and $x$ species, and as a result $f_{ee}\sim f_{xx}\sim f^{(FT)}$ and
similarly for the anti-neutrino flux factors. We leave a more detailed
analysis of possible closures to future work (Kneller et al., 2023).

###  4.2 EMU

EMU (Richers et al., 2021b) is a three-dimensional particle-in-cell neutrino
flavor transformation code that evolves Eqs. (5)-(6) individually for a large
number of computational particles. To evaluate the self-interaction part of
the Hamiltonian [Eq. (11)] we collect the contributions of each particle to
the background angular moments of the distribution using a second order shape
function, and interpolate the Hamiltonian from the grid to each particle using
the same second-order shape function. The advection terms are accounted for by
simply translating the position of each computational particle. The flavor
density matrix and positions of each computational particle are evolved with a
global fourth-order Runge-Kutta method. The snippets of EMU code that depend
on the number of neutrino flavors are automatically generated using sympy
(Meurer et al., 2017) to carry out symbolic matrix operations, simplify the
expressions, and output C++ code. This allows us to run simulations assuming
either two or three neutrino flavors. EMU is publicly available at Willcox &
Richers (2021).

Whereas FLASH is a moment method and only evolves two angular moments for each
flavor component of the neutrino distribution, EMU simulates particles moving
in many individual directions. The EMU results we present in this paper were
computed to 378 particles per cell corresponding to an angular resolution of
roughly 11 degrees, following the resolution tests in Richers et al. (2021a).

###  4.3 Initial and Boundary Conditions

We assign the flavor-diagonal values to the first two moments at every point
in the domain for a FLASH calculation. There are eight neutrino species with
four values (one energy density and three flux components), for a total of 32
initial values per cell which we need to assign. However, in practice we
always begin with identical moments for $\varrho_{xx}$ and
$\overline{\varrho}_{xx}$, since all heavy lepton neutrino and antineutrino
species are gathered in a single species $\nu_{\mathrm{other}}$ in Foucart et
al. (2016a) due to their very similar evolution in the absence of flavor
transformation.

Our calculations of the FFI will only include the self-interacting term for
the Hamiltonian-like operator in Eqs. (24), (25), and the anti-neutrino
counterparts. We set the vacuum and matter potentials to zero so as to focus
on the FFI, leaving the interesting physics cases of slow collective modes and
matter-neutrino resonances to future work. In an actual astrophysical object,
such as a CCSN or NSM, the vacuum potential would act to seed the flavor off-
diagonal elements as a function of path length and neutrino energy. Since we
simulate a local volume within a larger global system and thus have no
information about the advection of perturbed neutrinos into our domain, we
choose to take precise manual control over the initial seeds and exclude the
vacuum potential. We seed the off-diagonal flavor components with a
perturbation of $\mathcal{O}(10^{-6})$ compared to the diagonal components.
The scale $10^{-6}$ is chosen such that the growth in the off-diagonal
components begins in the linear regime. We have verified that starting with
even smaller perturbations does not change the outcome for FLASH calculations.
This is expected, as the growth should be the same in the linear regime, and
thus smaller initial perturbations only take longer and use more computing
resources. For the pattern of the perturbations, we use random numbers in each
cell in order to remain agnostic to scale of the initial perturbations that
would be present in nature.

Specifically, we use the following to seed the initial perturbations in the
off-diagonal components of the energy densities in FLASH

$\delta
E_{ab}(\mathbf{x})=10^{-6}\,p\,\max_{c}\\{N_{cc}\\}[A_{ab}(\mathbf{x})+iB_{ab}(\mathbf{x})],$
(47)

where $-1<A,B<1$ are uniform random numbers at each location $\mathbf{x}$,
$N_{cc}$ are the initial number density moments for the diagonal components,
and the Hermitian perturbation is only applied for $a\neq b$. For the flux
moment, we copy the energy density moment perturbation into the flux moment
and use flux factor vectors to weight the direction

$\delta\mathbf{F}_{ab}(\mathbf{x})=\delta
E_{ab}(\mathbf{x})\frac{\Sigma_{c}N_{cc}\mathbf{f}_{cc}}{\Sigma_{c}N_{cc}},$
(48)

implying the initial perturbations for the off-diagonal components of $E$ and
$\mathbf{F}$ are correlated. Analogous expressions exists for the anti-
neutrino moments.

In the EMU calculations, each particle is assigned a 4-momentum vector,
weight, and density matrix. The 4-momentum vectors are distributed uniformly
in space, but assigned initial weights and density matrices to approximate the
maximum entropy distribution [Eq. (34)] separately for each flavor. In this
way, the zeroth and first moments are reproduced under an appropriate angular
integration for each flavor-diagonal element of the density matrix. We impose
a random perturbation to the flavor off-diagonal elements of the density
matrix at the level of $10^{-6}$ and adjust the diagonal values accordingly to
preserve the length of each polarization vector.

The random numbers are determined at run time, so although the bulk properties
of the instability are expected to be related by the similar initial
conditions, the exact values in a particular cell have no correspondence
between EMU and FLASH calculations. We will make all comparisons in the
aggregate between the two sets of calculations. Furthermore, if we calculate
the energy density for EMU using Eq. (12) (where the angular integration
becomes a sum over particle index), we would expect the incoherent sum for
$\delta E_{ab}$ to be reduced by $\sqrt{n}$ for $n$ particles per cell,
implying an effectively smaller perturbation on the initial moments. To
reiterate, this only impacts the length of the linear phase of the instability
– not the growth rate or saturation properties.

In both FLASH and EMU calculations, we use a 3D cubic box with Cartesian
coordinates. The domain sizes and resolutions of the FLASH production
simulations are listed in Tables 1 and 3. We choose the domain size and cell
count so that we have the resolution to resolve the fastest growing mode in
the FFI, along with enough of a spatial domain to contain a few wavelengths of
that fastest growing mode. We do not know the properties of the fastest
growing mode a priori, so we perform convergence checks inline with the
presentation of the results. The simulation durations are generally longer
than the light crossing time of the domain, implying the initial
particles/densities will have advected out of the domain before the end of the
simulation. We implement periodic boundary conditions for both sets of
calculations implicitly assuming that the initial distribution is reasonably
approximated as periodic on scales larger than the simulation domain. We also
verify that changing the domain size does not impact the results.

##  5 Results

Our results comparing the ability of the two-moment method to reproduce the
fast-flavor instability are split into two parts. First, we consider in
Section 5.1 the three test problems in 3-dimensions that were previously
studied with EMU in Richers et al. (2021a). In Section 5.2 we move to
consideration of conditions extracted from a dynamical neutron star merger
simulation. We use the symbol ${\rm Im}(\Omega)$ to denote the growth rate of
$|N_{ex}|$ during instability. In addition, we use the symbol $|k|_{\rm max}$
to denote the fastest growing mode in the discrete fourier transform of
$N_{ex}$ during instability.

###  5.1 3D Test Problems

Name | $N_{ee}$ | $\overline{N}_{ee}$ | $\Sigma N_{(x)}$ | $\mathbf{f}_{ee}$ | $\overline{\mathbf{f}}_{ee}$ | $\mathbf{f}_{(x)}$ | $L$ | $N_{gp}$
---|---|---|---|---|---|---|---|---
$(10^{32}\,\mathrm{cm}^{-3})$ | $(10^{32}\,\mathrm{cm}^{-3})$ | $(10^{32}\,\mathrm{cm}^{-3})$ |  |  |  | (cm) |
Fiducial | $4.89$ | $4.89$ | 0 | $(0,0,1/3)$ | $(0,\phantom{/}0\phantom{3},-1/3)$ | $(0,0,0)$ | 8 | $128^{3}$
90Degree | $4.89$ | $4.89$ | 0 | $(0,0,1/3)$ | $(0,1/3,\phantom{-/}0\phantom{3})$ | $(0,0,0)$ | 8 | $128^{3}$
TwoThirds | $4.89$ | $3.26$ | 0 | $(0,0,\phantom{/}0\phantom{3})$ | $(0,\phantom{/}0\phantom{3},-1/3)$ | $(0,0,0)$ | 32 | $128^{3}$
Table 1: List of simulation parameters and initial conditions for the three 3D
test simulations. The first three columns show the number densities of each
anti/neutrino flavor. For clarity, the third column shows the sum of all four
heavy lepton anti/neutrino densities. Three-flavor simulations assume
$N_{\mu\mu}=\overline{N}_{\mu\mu}=N_{\tau\tau}=\overline{N}_{\tau\tau}=\Sigma
N_{(x)}/4$, while two-flavor simulations assume
$N_{xx}=\overline{N}_{xx}=\Sigma N_{(x)}/4$, where we assume that the other
half of the heavy-lepton neutrinos do not participate in flavor mixing. The
next three columns show the flux factor vectors, the norm of which are the
flux factors. The seventh column shows the length of each side of the domain
and the eighth column the number of grid points for the baseline simulation.
Figure 1: [Top] Polar representations of angular distributions for electron
neutrino (blue) and electron anti-neutrinos (red) for the three tests at the
beginning of the simulation. Blue (red) vectors indicate the net flux
direction. Purple vector is the sum of blue and red vectors. [Bottom] Curves
as given by Eq. (34) for $\nu_{e}$, $\overline{\nu}_{e}$ and the difference
(purple) as a function of polar angle $\theta$. Angular distributions for
$\nu_{x}$ and $\overline{\nu}_{x}$ are zero at the beginning of the
simulation. Lepton number crossings occur when purple line crosses $0$.

The three 3D test problems we consider are named as Fiducial, 90Degree, and
TwoThirds, all of which are described in detail in Richers et al. (2021b, a).
None of the three tests have analytic solutions222The Fiducial calculations of
Richers et al. (2021b) assume a slightly different angular distribution., so
the comparison is based on how well the moment method of FLASH can reproduce
the PIC results.

#### Test parameters

Table 1 gives the initial conditions for simulation parameters of the three
tests. The first three columns of Table 1 give the initial values of the
flavor-diagonal number density moment. All three tests start with non-zero
numbers of electron neutrinos and anti-neutrinos, and zero other-flavor
neutrinos. The fourth through sixth columns give the flux factor vectors.
Although these particular flux factor vectors need at most 2 dimensions to be
fully described, we stress that the calculations are three dimensional and
individual cells will generally develop flux moments where all three spatial
components are non-zero. The seventh column gives the side length of the
domain, and the eighth column the number of cells. Under the MEC, the initial
angular distributions of all three tests in Table 1 exhibit an ELN crossing
and are therefore unstable to FFC.

To visualize the geometry of these three tests, Fig. 1 shows the neutrino
angular distributions [Eq. (34)] for the electron neutrinos (blue) and anti-
neutrinos (red). As emphasized above, the MEC distributions are 3D, so the
polar representations are 2D cross-sectional slices. We take the slices such
that the maximum values of the distributions lie in the same plane. The polar
plots in the upper panels show a more intuitive representation of the
magnitude of the distribution in different directions, but the size and depth
of the ELN crossings are more apparent in the standard plots on the lower
panel. Indeed, the vectors originating from the origin on the polar plots show
the peak direction of the distributions, with a reversed direction for anti-
neutrinos (in order to give a sense of the net flux of lepton number). The sum
of the blue and red vectors, shown in purple, is coincident with the red and
blue ones for the Fiducial test, and coincident with the red one for the
TwoThirds test. The 90Degree test shows vividly the vector sum. For all three
tests, we orient the polar plane so that the purple vector points in the
rightwards direction, and indicate this direction as $\theta=0$ in the lower
plots. As in the polar plots, the blue and red curves give the electron
neutrino and the negative value of the electron anti-neutrino distributions.
Here, though, the purple curve is the sum of the blue and the red curves, and
gives the ELN distribution, properly normalized by the sum of the energy
density moments. As can clearly be seen in all three tests, the purple curves
cross the horizontal axis implying a lepton number crossing.

Figure 2: Density matrix elements versus time for the three tests. Two-flavor
FLASH quantities are plotted in red. Two- (Three-) flavor EMU quantities are
plotted in solid (dashed) black. All quantities are averaged over the spatial
domain, and in addition over particle number for EMU. [Top] Plotted is the
$ee$ component of the number density moment flavor-matrix $N$ (i.e., number
density of $\nu_{e}$) scaled by $N_{ee}(t=0)$. [Bottom] Plotted is the
magnitude of the off-diagonal component of $N$ scaled by $N_{ee}(0)$. For
3-flavor EMU calculations, we take the $e\mu$ component of $N$.

#### Time evolution and FFI

We show the time-evolution of the domain-averaged values of
$N_{ee}(t)/N_{ee}(0)$ (top) and $|N_{ex}(t)|/N_{ee}(0)$ (bottom) in Fig. 2 for
FLASH, 2-flavor EMU, and 3-flavor EMU simulations. The saturation times
$t_{\mathrm{sat}}$ are at the peak of the $|N_{ex}|/N_{ee}(0)$ curves, which
we use to line up results from simulations with different initial
perturbations. In addition, we plot a horizontal green line on the top panel
to indicate the expected number of electron neutrinos in a _2-flavor_
calculation if the system were to completely mix flavor.

In all three tests, the FLASH simulations exhibit fast flavor instability with
a growth rate very similar to the true value. Considering the red (FLASH)
curves in Fig. 2, there exists a period of exponential growth in
$\langle|N_{ex}|\rangle$, evidencing one of the defining characteristics of
the FFI. $\langle|N_{ex}|\rangle$ continues to grow until the off-diagonal
magnitude reaches the same order of magnitude as the initial electron-flavor
number density moment. When $|N_{ex}|\lesssim N_{ee}$ at saturation, there are
rapid oscillations in the diagonal components, evidencing the other defining
characteristic of FFI. Saturation is a nearly instantaneous event with
decoherence, i.e., decreasing $|N_{ex}|$, succeeding the rapid oscillations.
The decoherence continues as oscillations damp, with an end result of $\langle
N_{ee}\rangle$ approaching an asymptote at a value less than the starting
condition. In summary, the results presented in Fig. 2 are quite remarkable in
that even though instability criterion in FFI depends upon angular crossings
of the ELN, the two-moment method accurately showcases the growth of the FFI
without access to crossing information. Of course, crossings are implied by
the MEC distribution used to generate the closure relation, but the MEC
distribution is nowhere explicitly used in the code.

The growth rate in the FLASH simulations is quantitatively very similar to
that in the full EMU simulations, and even the final asymptotic neutrino
distributions match well in certain cases. The Fiducial and 90Degree tests
show strong agreement between the two methods. We see nearly identical growth
rates for both tests, with FLASH producing a slightly higher ${\rm
Im}(\Omega)$. Specifically for the Fiducial test case, ${\rm
Im}(\Omega)=7.08\times 10^{10}\,{\rm s}^{-1}$ in FLASH, as compared to ${\rm
Im}(\Omega)=6.29\times 10^{10}\,{\rm s}^{-1}$ in EMU, for a difference of
$\sim 10\%$. Results are similar for the 90Degree test, with ${\rm
Im}(\Omega)=5.38\,(4.36)\times 10^{10}\,{\rm s}^{-1}$ in FLASH (EMU). Also,
the asymptotic values for $\langle N_{ee}\rangle$ are nearly the same, with
differences of $\sim 1\%$ for both tests. We nevertheless see differences
between these two tests in the saturation and decoherence periods. The
oscillations in the top panels of Fig. 2 for EMU appear to have larger
amplitudes and persist longer than those of FLASH.

The small amplitude oscillations for FLASH are especially evident in the
TwoThirds test case. Here, we see a noticeable difference between FLASH and
EMU for the asymptotic values of $\langle N_{ee}\rangle$. The growth rate is
faster for FLASH by $\sim 40\%$ and the loss of coherence falls off faster.
There is a smaller amount of time when $|N_{ex}|\lesssim N_{ee}$, and thus
less oscillations in the flavor-diagonal term. The result is an asymptotic
value which is $\gtrsim 10\%$ of $N_{ee}(0)$.

For the TwoThirds test, we speculate that the reason the moment calculations
do not asymptote at large times to the same value of $\langle N_{ee}\rangle$
as found by EMU is due to our imposition of the MEC. Recall that for the FLASH
calculations, we use the quantum implementation of the MEC at every time step
and substep of the evolution. In contrast, the EMU calculations use Eq. (34)
only when generating the initial conditions, and the future evolution depends
directly on the general distribution. There is no guarantee that the neutrino
distributions in EMU follow the MEC at any point except for initialization.

Although we have argued above that the use of the MEC in FLASH necessarily
restricts the shape the distributions may take during flavor evolution, there
does exist the striking convergence between FLASH and EMU of $\langle
N_{ee}\rangle$ in the asymptotic limit for the Fiducial and 90Degree tests.
This is not a coincidence, but rather a result of the symmetry of both of
these tests. Initially, the system contains both $CP$ and rotational
symmetries. The MEC is agnostic to $CP$ but does preserve the rotational
invariance for constant flavor-traced flux factors. As the energy density
moment is equal for $E_{ee}$ and $\overline{E}_{ee}$, and the initial neutrino
distributions are rotations of the anti-neutrino ones: an ELN crossing is
inevitable. The initial conditions and conservation of 3-momentum ensures that
neutrinos and anti-neutrinos will never have the same flux factor vectors at
any point in the test calculations. As our system of equations is $CP$
symmetric (except for the initial conditions in the flux factors), we expect
any flavor transformation for $E_{ee}$ to be accompanied by a commensurate
change in $\overline{E}_{ee}$. Because 3-momentum is conserved, the flux
factors are invariant and the ELN crossing persists to all times. We have
numerically verified that indeed $\overline{E}_{ee}$ mirrors the evolution of
$E_{ee}$ and an ELN crossing exists in perpetuity. In other words, the
distributions shown in the top panels of Fig. 1 only scale in radial
coordinate during their evolution. However, the results in Fig. 2 clearly show
a stable system post saturation. For either the Fiducial or 90Degree system to
become stable, an XLN crossing must develop, cancelling the omnipresent ELN
one. Furthermore, the $\nu_{x}$ and $\overline{\nu}_{x}$ distributions have
the same vector flux factors and use the same flavor-traced flux factor,
implying those distributions are identical to the ones in the top panels of
Fig. 1 except for a difference in the radial coordinate. In the presence of
non-trivial ELN and XLN crossings, a zero net lepton number at all angles
requires identical energy, flux, and pressure moments for the $xx$ components
as compared to their $ee$ counterparts – implying near flavor equilibration.
Even if the distributions do not follow Eq. (34) and the MEC, the symmetry of
the system guarantees that $\langle N_{ee}\rangle$ must converge to $50\%$ of
the flavor trace [equivalent to $N_{ee}(0)/2$] in both the Fiducial and
90Degree tests. This need not be the situation in the TwoThirds test case as
the system neither exhibits $CP$ or rotational symmetry. Here, the MEC is not
an accurate representation of the distributions at later times, and as a
result, the FLASH and EMU calculations show a stark divergence. Discrepancies
between moment and multi-angle methods were also seen in Myers et al. (2022).

Figure 3: Magnitude of the discrete Fourier transform of $N_{ex}$ for all
three tests plotted against wavenumber. The DFTs are calculated at a time
prior to saturation. The DFTs for the Fiducial and 90degree test cases are
taken $0.1\,{\rm ns}$ before saturation, and TwoThirds $0.35\,{\rm ns}$ before
saturation. 2-flavor (3-flavor) EMU simulations correspond to the black solid
(dashed) curve; FLASH simulations to the red solid curve. The two light-red
curves correspond to different resolution tests for FLASH. The light-red solid
curve has half the box-side length and half the number of grid points compared
to columns 7 and 8 in Table 1. The light-red dashed curve has the same box-
side length and half the number of grid points.

#### Fourier space analysis

We have discussed averages of the number density moment when presenting Fig.
2. In Fig. 3 we show information on the structure in the simulation domain by
using Discrete Fourier Transforms (DFTs). The horizontal axes give the
wavenumber $k$, and the vertical axis the magnitude of the DFT of the complex
flavor-off-diagonal number density moment $N_{ex}$, normalized by the flavor
trace. We will refer to wavenumber values as “modes”. The solid red, solid
black, and dashed black lines all correspond to the same simulations as Fig.
2. In lighter shades of red we have plotted DFTs from two additional FLASH
calculations of the same test cases. The light-red solid curve is from a
simulation with the same number of grid points per cm but with a box-side-
length of half the original simulation compared to the values in columns 7 and
8 of Table 1; the light-red dashed curve is from a calculation with the same
domain size, but half the number of grid points per cm and a smaller maximum
value of $k$. The three DFTs in each panel of Fig. 3 are all from a time
before saturation during the growth period: $\sim 0.1\,{\rm ns}$ before
saturation for the Fiducial and 90Degree cases; $\sim 0.35\,{\rm ns}$ before
saturation for the TwoThirds case. While similar, the times of the snapshots
used in the DFTs are not exactly equal between different calculations so the
values of $\widetilde{N}_{ex}$ cannot be compared across either the
simulations or the resolution tests. For this reason, comparisons should be
restricted to within an individual calculation, i.e., the relative heights of
peaks.

The DFTs show the scales, via wavenumber $k$, where there exists a sinusoidal
pattern in the flavor off-diagonal number density moment. This superposition
of sinusoids need not have growing amplitudes for each mode. A priori, only
one mode is necessary to explain the growth phase in Fig. 2. However, during
the growth phase, all modes in Fig. 3 do indeed grow in power until
saturation, implying there are many unstable modes in the system.

All three tests show a discernible peak in the dark red curves of Fig. 3. Soon
after the simulations begin the DFT exhibits a peak with an associated wave
number as evidenced in Fig. 3. The peak remains at that location in $k$,
although with growing height, until saturation. The DFTs for the resolution
tests show similar behavior in the peak position and growth phase, indicating
that the dark red curve for the simulation is indeed spatially resolved. We
call the wavenumber at this peak the fastest growing mode $|k|_{\rm max}$. The
growth rate in Fig. 2 is linked to the fastest growing mode via a dispersion
relation, with details provided in Froustey et al. (2023).

The FLASH and EMU calculations both have discernible peaks with similar
fastest growing modes. The wave number of the fastest growing modes for FLASH
are slightly larger, reflecting a smaller scale. For example, $|k|_{\rm
max}=3.93\,(3.14)\,{\rm cm}^{-1}$ for FLASH (EMU) in the Fiducial test case,
and $|k|_{\rm max}=3.14\,(2.36)\,\mathrm{cm^{-1}}$ in the 90Degree case. Also,
it appears that the noise floor of the DFT is larger for FLASH, or
equivalently, there exists relatively less power in the fastest growing mode.
Lastly, there are a few harmonics visible in FLASH but not present in EMU.
This is true for all three test cases, and more pronounced for the TwoThirds
case. These harmonics, however, only crest slightly above the noise floor.

Name | $\langle|N_{ex}|/N_{ee}(0)\rangle|_{t=t_{\rm sat}}$ | $\langle N_{ee}/N_{ee}(0)\rangle|_{t\rightarrow\infty}$ | ${\rm Im}(\Omega)$ | $|k|_{\rm max}$
---|---|---|---|---
$(10^{10}\,\mathrm{s}^{-1})$ | $(\mathrm{cm}^{-1})$
Fiducial |  |  |  |
FLASH | $0.314$ | $0.517$ | $7.08$ | $3.93$
EMU (2f) | $0.333$ | $0.506$ | $6.29$ | $3.14$
90Degree |  |  |  |
FLASH | $0.281$ | $0.510$ | $5.38$ | $3.14$
EMU (2f) | $0.303$ | $0.516$ | $4.36$ | $2.36$
TwoThirds |  |  |  |
FLASH | $0.248$ | $0.627$ | $1.95$ | $1.77$
EMU (2f) | $0.214$ | $0.771$ | $1.20$ | $1.37$
Table 2: Numerical results for FLASH and EMU (2f) calculations for the three
3D test simulations. First column gives the ratio $\langle\lvert
N_{ex}\rvert/N_{ee}(0)\rangle$ when $t=t_{\rm sat}$. Second column gives the
asymptotic ratio $\langle N_{ee}/N_{ee}(0)\rangle$ post-saturation. Third
column gives the growth rate ${\rm Im}(\Omega)$ when the system is unstable in
units of $10^{10}\,{\rm s}^{-1}$. Last column gives the fastest growing mode
in the domain $|k|_{\rm max}$ in units of ${\rm cm}^{-1}$.

In summary, Table 2 gives numerical results of FFC to compare between FLASH
and the 2-flavor EMU calculations for all three tests. The values in columns
one through three are deduced from Fig. 2 and are the following, respectively:
the maximum value in the bottom panels of Fig. 2; the asymptote in the top
panel of Fig. 2; and the slope of the line (in semi-log space) in the bottom
panels of Fig. 2. The fourth column gives the value of $k$ at the peak of the
DFT in Fig. 3. All tests show the FLASH calculations have larger values of
${\rm Im}(\Omega)$ and $|k|_{\rm max}$ compared to EMU. However, even with the
different growth rates, oscillations begin once the average value of
$|N_{ex}|$ reaches quite similar values for each test. This value, in the
first column, is similar within a given test but not uniform across all three
tests.

###  5.2 Neutron Star Merger

Our next set of simulations use initial conditions extracted from the three-
dimensional neutron star merger simulation of Foucart et al. (2016a). This
simulation is general relativistic, but simulating neutrino oscillations in
curved spacetime is beyond the scope of this work. Appendix B describes our
procedure on transforming the distributions defined in a general spacetime
metric to distributions defined in an orthonormal tetrad comoving with the
fluid. In this frame, the construction of a flavor transformation simulation
is more intuitive, since we can treat the spacetime as locally flat.

Figure 4: NSM crossing information from Foucart et al. (2016a). The snapshot
is taken $5\,{\rm ms}$ post merger. Red pixels indicate locations where an ELN
crossing exists. The three symbols (black cross, blue star, green circle)
indicate the locations for the three flavor-transformation simulations we
consider in Sec. 5.2. From outside to inside, the green contours indicate
matter densities of {$10^{11},10^{12},10^{13},10^{14}$} $\mathrm{g.cm}^{-3}$.

We analyze and simulate neutrino distributions at a selection of points in the
polar slice of a snapshot at $5\,{\rm ms}$ post merger shown in Fig. 4. Green
contours give matter densities of {$10^{11},10^{12},10^{13},10^{14}$}
$\mathrm{g\,cm}^{-3}$ and the inner contours show the position of the central
compact object at the center of the domain. The red pixels indicate where an
ELN crossing exists according to Eq. (41). White pixels indicate that no such
ELN crossing exists, although these regions are still subject to flavor
transformation via other processes (e.g., the matter-neutrino resonance) or
advection of flavor-transformed distributions into those regions of space. We
select three points to simulate in FLASH and EMU, indicated by the black
cross, blue star, and green circle, to model regions above the accretion disk,
within the disk, and just outside of the compact object, respectively. The
black cross is the same point detailed in Grohs et al. (2022).

Name | $N_{ee}$ | $\overline{N}_{ee}$ | $\Sigma N_{(x)}$ | $\mathbf{f}_{ee}$ | $\overline{\mathbf{f}}_{ee}$ | $\mathbf{f}_{(x)}$ | $L$ | $N_{gp}$
---|---|---|---|---|---|---|---|---
$(10^{32}\mathrm{cm}^{-3})$ | $(10^{32}\mathrm{cm}^{-3})$ | $(10^{32}\mathrm{cm}^{-3})$ |  |  |  | (cm) |
NSM 1 | 14.22 | 19.15 | 19.65 | $\begin{pmatrix}\phantom{-}0.0974\\\ \phantom{-}0.0421\\\ -0.1343\end{pmatrix}$ | $\begin{pmatrix}\phantom{-}0.0723\\\ \phantom{-}0.0313\\\ -0.3446\end{pmatrix}$ | $\begin{pmatrix}-0.0216\\\ \phantom{-}0.0743\\\ -0.5354\end{pmatrix}$ | 7.87 | $128^{3}$
NSM 2 | 23.29 | 28.53 | 60.11 | $\begin{pmatrix}\phantom{-}0.0086\\\ -0.0174\\\ -0.1635\end{pmatrix}$ | $\begin{pmatrix}\phantom{-}0.0070\\\ -0.0142\\\ -0.2338\end{pmatrix}$ | $\begin{pmatrix}-0.0476\\\ -0.0231\\\ -0.2679\end{pmatrix}$ | 8.27 | $256^{3}$
NSM 3 | 28.80 | 37.42 | 19.32 | $\begin{pmatrix}\phantom{-}0.0004\\\ -0.0033\\\ \phantom{-}0.0044\end{pmatrix}$ | $\begin{pmatrix}\phantom{-}0.0003\\\ -0.0025\\\ -0.1306\end{pmatrix}$ | $\begin{pmatrix}-0.0008\\\ -0.0051\\\ -0.1292\end{pmatrix}$ | 5.80 | $512^{3}$
Table 3: List of baseline simulation parameters for the FLASH NSM simulations.
Column labels are the same as Table 1. Note that all corresponding EMU
simulations were run with the same parameters, but using $N_{gp}=128^{3}$ grid
cells due to the longer wavelength of the fastest growing mode.  Figure 5:
Polar and cartesian representations of the initial $\nu_{e}$ and
$\overline{\nu}_{e}$ distributions for the NSM points. Plotting conventions
are the same as Fig. 1.

Table 3 and Fig. 5 are the NSM analogs to Table 1 and Fig. 1 of Sec. 5.1, and
use the same notation and plotting conventions. Note that the
orthogonalization procedure in App. B is location-dependent and as a result
the directions of the fluences in Table 3 cannot be compared to one another
between points. In other words, for this particular study our flavor
transformation simulations are restricted to the local area of each point and
do not affect one another through advection. However, the neutrinos for points
1 and 2 are generally moving upward, so the leftward direction in the polar
plots in Fig. 5 roughly correspond to the $\hat{z}$ direction in Fig. 4. The
leftward direction in the polar plot for point 3 in Fig. 5 roughly corresponds
to the $\hat{x}$ direction in Fig. 4. Note that while there are healthy ELN
crossings in points 1 and 3, the crossings in point 2 are quite tenuous, as
would be expected given the very thin band of instability just above the
compact object in Fig. 4.

Figure 6: Domain averaged components for the number density moment plotted
against time measured from the point of saturation. Panels and curve
conventions are similar to Fig. 2 and simulation computational parameters are
given in the last two columns of Table 3. Gray and light red lines give
results from resolution tests and are dependent on the NSM point. For NSM 1,
the solid (dashed) gray lines are for EMU 2-flavor (3-flavor) simulations with
half the side-length and half the number of grid points per side, for an
identical grid spacing. The solid medium-red line also is for a FLASH
simulation with half the side-length and half the number of grid points. In
addition, the dashed light-red line is for a simulation with half the side-
length but the same number of grid points, for half the grid spacing. For the
3-flavor EMU calculation in NSM 1, $ex=e\mu$. NSM 2 and 3 follow identical
conventions for resolution testing compared to one another. Gray lines are for
2-flavor EMU simulations with half the domain size and half the number of grid
points per side. Medium-red FLASH simulations are also half the domain size
and half the number of grid points per side. The light-red FLASH calculation
is for the same domain side-length, but only half the number of grid points
resulting in twice the grid spacing of the baseline simulation.

#### Time evolution and FFI

Figure 6 shows the time-evolution of the neutrino number density moment for
all three NSM points. Contrary to the previous test cases, the conditions in
the NSM dictate non-zero initial distributions of heavy lepton neutrinos. This
results in different 2-flavor complete mixing lines for each simulation, shown
in green in Fig. 6. The lighter-opacity lines are for different resolution
tests (see descriptions in caption). For illustrative purposes, we also
include a three flavor EMU simulation for the first NSM point and plot the
$e\mu$ component in the bottom panel.

In all three NSM points, we see growth, saturation, and decoherence phases as
we did in Sec. 5.1 and Fig. 2. Growth of $|N_{ex}|$ begins soon after the
start of the simulation and proceeds until $\langle|N_{ex}|\rangle\sim
0.1\langle N_{ee}\rangle$ in both FLASH and EMU calculations. Rapid
oscillations develop and effect a decrease in $\langle N_{ee}\rangle$ towards
complete flavor-mixing. Notice that for the first NSM point treated by FLASH,
the conditions are such that $\langle N_{ee}\rangle$ falls below the complete-
mixing green line and asymptotes to a value less than $50\%$ of the flavor
trace. This is the case for all three resolution tests, including the light-
red dashed curve with half the gird spacing compared to the baseline
simulation. The EMU results also briefly dip below the 50% line, and it seems
that the more rapid decoherence in the moment method halts the flavor
transformation before it can oscillate back up. In all other calculations
(FLASH and EMU) $\langle N_{ee}\rangle$ remains above the green line at all
times. Decoherence enters after saturation in much the same manner as the
three test calculations in Sec. 5.1. In the first and second NSM points, the
baseline and resolution tests for FLASH begin to lose convergence in the
decoherence phase. The divergence occurs well after saturation and at a point
where $\langle N_{ee}\rangle$ has reached a steady-state value. The FLASH
baseline and resolution tests for the third NSM point maintain convergence
longer – a result of this set of calculations having smaller grid spacings
compared to the other two points.

We identify some general trends in the FLASH and EMU results. FLASH generally
shows faster growth rates, faster decoherence fall-offs, and less oscillations
in the $N_{ee}$ moment, similar to the test cases in the previous section. The
discrepancies are particularly apparent for the second NSM point. This point
is unique in that the distribution described by the MEC is only marginally
unstable. This type of condition is expected to lead to slower growth rates,
less total flavor transformation, and more dependence on details in the small
angular region between the ELN crossing points (e.g., Richers et al. 2021b;
Bhattacharyya & Dasgupta 2022). Specifically for FLASH (EMU), ${\rm
Im}(\Omega)=5.19\,(1.10)\times 10^{10}\,{\rm s}^{-1}$ for this point. The
bottom panels of Fig. 5 shows that ELN crossings are initially present for all
three points, but the crossing is most shallow for the second point. The two
moment method plus MEC we are employing in FLASH is not able to capture the
FFI behavior as well for this scenario as it is for the more pronounced ELN of
the first point, and seems to behave as if there was a more significant
instability than present in the detailed angular distribution, although we
again emphasize that FLASH still demonstrates a characteristic evolution
pattern for a fast-flavor-unstable distribution in general.

Figure 7: Volume rendering of contours of the phase of $N_{ex}$, $\phi_{ex}$,
for NSM 1 (black cross in Fig. 4). Blue, white, and red contours correspond to
$\phi_{ex}=-\pi/2,0,\pi/2$, respectively. The four panels are at four
different times and roughly correspond to different phases of FFC. From left
to right: $t=0.01\,{\rm ns}$ (initial conditions); $t=0.18\,{\rm ns}$ (growth
phase); $t=0.21\,{\rm ns}$ (saturation point); and $t=0.58\,{\rm ns}$
(decoherence phase).

We show in Fig. 7 3D-volume-renderings of the FLASH simulation for the first
NSM point at four different simulation times: $t=0.01\,{\rm ns}$;
$t=0.18\,{\rm ns}$; $t=0.21\,{\rm ns}$; and $t=0.58\,{\rm ns}$. The contours
in each panel are for the phase, $\phi_{ex}$, of the complex number $N_{ex}$.
The spatial structure in $\phi_{ex}$ reflects the phase of the growing mode,
and so reflects the three-dimensional structure of the peak of the DFT during
the linear phase and a combination of the persisting mode structure and random
decoherence after the instability saturates. We plot three contours for the
phase: $\phi_{ex}=-\pi/2$ (blue); $\phi_{ex}=0$ (white); and $\phi_{ex}=\pi/2$
(red). The first panel shows a time close to the start of the simulation,
where little flavor-transformation has occurred and the contours are close to
the initial conditions of the random distributions with no structure [Eq.
(47)]. We can see some structure in the second panel during the growth phase,
where the distance between planar structures reflects the wavelength of the
fastest growing mode and the planes are roughly perpendicular to the direction
of the net ELN flux. The phases are distorted when the evolution becomes
nonlinear as the instability saturates in the third panel. The last panel is
during the decoherence phase when the flavor field is no longer unstable, yet
there still exists structure in $\phi_{ex}$. The pattern seen here in Fig. 7
is qualitatively similar to that seen in Fig. (2) of Richers et al. (2021a),
albeit for a different simulation that employed the PIC method. Nevertheless,
the similarity in the growth of structure of $\phi_{ex}$ again shows that the
moment method reproduces many features of the FFI on large and small scales.

Figure 8: The magnitude of the discrete Fourier transform of $N_{ex}$ in the
three NSM simulations plotted against wavenumber. DFTs are at a time prior to
saturation. For NSM 1, both FLASH and EMU DFTs are $0.1\,{\rm ns}$ before
saturation. For NSM 2, the DFT for FLASH occurs $0.5\,{\rm ns}$ before
saturation and $0.1\,{\rm ns}$ for EMU. For NSM 3, the DFT for FLASH occurs
$0.15\,{\rm ns}$ before saturation and $0.05\,{\rm ns}$ for EMU.

#### Fourier space analysis

Figure 8 gives the DFTs for the three NSM points at $0.1\,{\rm ns}$ before
saturation, roughly corresponding to the second panel in Fig. 7. Because the
growth rates are different for the second (third) points, we pick a time
before saturation of $0.5\,{\rm ns}$ ($0.15\,{\rm ns}$) for FLASH, and
$0.1\,{\rm ns}$ ($0.05\,{\rm ns}$) for EMU, all of which allow us to most
clearly capture the fastest growing mode before it begins to nonlinearly
couple to other modes. For the first NSM point, we find $|k|_{\rm
max}=6.39\,{\rm cm}^{-1}$, corresponding to a wavelength of $1.02\,{\rm cm}$
and matching the distance between planar structures in the second panel of
Fig. 7 (domain size $L\sim 8\,{\rm cm}$).

Similar to the results of Sec. 5.1, the DFTs from the FLASH simulations show
fastest growing modes with smaller wavelengths, higher noise floors, and
visible harmonics for all three simulation points. The EMU simulations were
run with a larger grid cell size to most optimally resolve the larger unstable
wavelengths, resulting in a DFT that is cut off at smaller maximum $k$.
Despite the differences between the FLASH and EMU calculations, we emphasize
that all of the FLASH simulations once again reflect characteristic FFI
behavior with discernible peaks with reasonable values of the fastest growing
mode.

Name | $\langle|N_{ex}|/N_{ee}(0)\rangle|_{t=t_{\rm sat}}$ | $\langle N_{ee}/N_{ee}(0)\rangle|_{t\rightarrow\infty}$ | ${\rm Im}(\Omega)$ | $|k|_{\rm max}$
---|---|---|---|---
$(10^{10}\,\mathrm{s}^{-1})$ | $(\mathrm{cm}^{-1})$
NSM 1 |  |  |  |
FLASH | $0.158$ | $0.643$ | $8.09$ | $6.39$
EMU (2f) | $0.178$ | $0.743$ | $5.58$ | $4.79$
NSM 2 |  |  |  |
FLASH | $0.0845$ | $0.837$ | $5.19$ | $6.08$
EMU (2f) | $0.0640$ | $0.960$ | $1.10$ | $3.80$
NSM 3 |  |  |  |
FLASH | $0.181$ | $0.609$ | $10.7$ | $13.0$
EMU (2f) | $0.170$ | $0.831$ | $4.15$ | $6.50$
Table 4: Results for FLASH and EMU (2f) calculations for the three 3D NSM
simulations. Column labels are the same as Table 2.

To summarize, Table 4 gives numerical results of FFC to compare between FLASH
and the 2-flavor EMU calculations for the three NSM simulations. The results
in Table 4 are presented in the same way as Table 2. Like the Fiducial,
90Degree, and TwoThirds tests, the FLASH calculations for the three NSM points
have larger values of ${\rm Im}(\Omega)$ and $|k|_{\rm max}$ compared to EMU,
but still exhibit reasonable behavior characteristic of the FFI.

##  6 Conclusions

Core collapse supernovae and merging neutron stars are complex systems that
require the melding of many different physics aspects including
magnetohydrodynamics, general relativity, equation of state physics, and
neutrino physics. When neutrino moment methods are used currently in large
scale simulations, they consist of classical neutrino physics and typically
employ two angular moments with a closure for each neutrino species. In this
work, we have extended the moment method framework in the context of the FLASH
code to take into account neutrino flavor transformation in a two-moment
scheme with a quantum closure. While at present we have tested the flavor
transformation alone, the success of the moment method in modeling many
features of the fast flavor instability lends confidence to incorporating it
into large-scale multi-physics simulations in the future.

We found that neutrino transformation behavior is well captured in a number of
test problems, including vacuum oscillations, Mikheyev-Smirnov-Wolfenstein
(MSW) resonance, and bipolar oscillations. Additionally, we approximately
reproduced the results of three multidimensional PIC simulations in the
literature at a fraction of the computational cost. For those tests, in the
key metric of the final electron neutrino number density, we find very good
agreement among the first two ($\sim 1\%$) and more qualitative agreement with
the last ($\sim 15\%$). Similarly, the two-moment simulations were able to
almost exactly match instability growth rates in the first two although our
moment method has somewhat faster growth in the last. The moment method also
shows fastest growing modes that peak at a similar, but slightly higher wave
number than in the PIC calculations. Further analysis is required to tease out
the details of the dispersion relation under the two-moment approximation in
regions of instability (Froustey et al., 2023) and compare with numerical
simulation.

We then performed multidimensional simulations of the FFI in three separate
neutrino angular distributions taken from a full-scale classical simulation of
an NSM. In key quantities, we found similar levels of agreement as we found in
the TwoThirds test. Specific differences include a lower electron neutrino
number density at saturation with the moment method than with EMU, and a
faster growth rate with the moment method. While we always found qualitative
agreement, the different methods naturally show the largest deviations when
the distribution is only marginally unstable. Quantitative agreement is best
with deep crossings and tends to worsen in the case of shallow crossings. This
shortcoming will be important to improve upon in future advances of the
algorithm, since the rapid onset of the FFI is likely to drive ELN crossings
to remain shallow in astrophysical environments.

Our two-moment algorithm is based on an extension of the classical MEC
relevant to quantum neutrino transport (Richers, 2020). Since the closure
accounts for missing physics from the unevolved higher-order moments, we
anticipate that future efforts to develop quantum closures will improve the
agreement between this two-moment method and more exact methods. To go beyond
the simple prescription used here requires using all of the components of $E$
and $\mathbf{F}$ in a basis-independent way as suggested by Kneller et al.
(2023).

Because of the small scales on which flavor transformation occurs, our method
is at present still too computationally expensive to directly place in a full-
scale CCSN or NSM simulation. However, we anticipate that in the future,
methods can be developed to incorporate the very small-scale physics of flavor
transformations into large-scale simulations. Nevertheless, by virtue of only
following two moments, this method is substantially computationally cheaper
than exact ones that evolve neutrino distributions along hundreds or thousands
of directions. Given the successes of capturing many features of the FFI in
this work, we believe that moment methods can complement the higher fidelity
methods such as PIC and multi-angle. For example, moment methods could be
useful in doing faster realistic calculations when searching configuration
space (Johns & Nagakura, 2021), with follow-up post-processing being done by
PIC or multi-angle codes. Machine learning is also a possibility, with moments
being used to train the algorithms (Abbar, 2023).

Despite the caveats above, we find that our angular-moment implementation of
the QKEs reasonably and effectively captures the complex and confounding
phenomenon of neutrino flavor transformation in conditions which are plausible
in the environments of CCSNe and binary NSMs. Specifically, in anisotropic
conditions where the angular neutrino distributions exhibit a lepton number
crossing, the two moments of an M1 transport scheme manifest the phases of
fast flavor conversion: exponential growth during unstable conditions; peak
saturation of the off-diagonal density matrix component; rapid flavor
oscillations of the diagonal components; and post-saturation decoherence with
subsequent freeze-out.

There are many possible extensions that one can make to the work presented
here. Improvements on the closure, the addition of the collision matrix, and
an extension to three flavors are a few. Finally, while we have included
advection in our simulations, the inclusion of both advection and flavor
mixing in the context of a large scale simulation will likely alter the
angular distributions of the neutrinos Padilla-Gay et al. (2021). These
improvements will widen the conditions for which this method can be used and
allow us to probe other predicted phenomena, such as collisional
instabilities, bipolar oscillations and matter neutrino resonance transitions.

Including flavor transformation in 3D general-relativistic-
magnetohydrodynamics astrophysics simulations is a major computational
challenge for multi-messenger-astrophysics theory. Our moment-method flavor
calculations offer a contribution to this important field of study.

##  7 Acknowledgments

We thank George Fuller and MacKenzie Warren for useful conversations. EG, JPK,
and GCM are supported by the Department of Energy Office of Nuclear Physics
award DE-FG02-02ER41216. EG acknowledges support by the National Science
Foundation grant No. PHY-1430152 (Joint Institute for Nuclear Astrophysics
Center for the Evolution of the Elements). SMC is supported by the U.S.
Department of Energy, Office of Science, Office of Nuclear Physics, Early
Career Research Program under Award Number DE-SC0015904. JF is supported by
the Network for Neutrinos, Nuclear Astrophysics and Symmetries (N3AS), through
the National Science Foundation Physics Frontier Center Grant No. PHY-2020275.
SR is supported by a National Science Foundation Astronomy and Astrophysics
Postdoctoral Fellowship under award AST-2001760. FF is supported by the
Department of Energy, Office of Science, Office of Nuclear Physics, under
contract number DE-AC02-05CH11231 and by the National Science Foundation
through grant AST-2107932. The calculations presented in this paper were
undertaken, in part, on the _Payne_ machine at North Carolina State University
which is supported in part by the Research Corporation for Science
Advancement. This research used resources of the National Energy Research
Scientific Computing Center (NERSC), a U.S. Department of Energy Office of
Science User Facility located at Lawrence Berkeley National Laboratory,
operated under Contract No. DE-AC02-05CH11231 using NERSC award ERCAP0021631.
This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Advanced Scientific Computing Research and Office
of Nuclear Physics, Scientific Discovery through Advanced Computing program
under Award Number DE-SC0017955. The authors would like to acknowledge the use
of the following software: Matplotlib (Hunter, 2007), Numpy (van der Walt et
al., 2011), and SciPy (Jones et al., 2001).

##  Appendix A 1D Test Problems

To test the algorithms for simultaneous flavor-transformation and advection we
implemented in FLASH, we have designed and conducted three tests in 1D to
quantify the numerical accuracy and precision of the code. This appendix
contains the results of those tests, which we term: Vacuum; MSW; and Bipolar
Oscillations. We compare to analytic solutions for the vacuum and MSW tests,
and to an analytic prediction for the period of the bipolar test.

###  A.1 Vacuum Oscillations

In this test problem we consider only the vacuum potential $H_{V}$ in the QKEs
of Eqs. (24) and (25), but with a transformation into spherical coordinates.
Neutrinos and anti-neutrinos are emitted isotropically from a point source at
the origin. The domain has an outer radius of $4\,{\rm km}$. We use a
spherical geometry with an inner radius of $r_{0}=100\,{\rm m}$ to avoid
divide by zero errors at the origin, and assume the flux factors at $r_{0}$
are all unity as the neutrinos stream outwards in the radial direction. The
number of grid points is 640, all evenly-spaced in radius. Initially, all
cells have the same density of the four particle species scaled by $1/r^{2}$.
The flux factors are set to 1.0. The energy densities of the off-diagonal
species are set to zero. We have set the outer boundary condition to an
outflow, and the inner boundary to that of an oven emitting
$\nu_{e},\overline{\nu}_{e},\nu_{x},$ and $\overline{\nu}_{x}$. At time $t=0$
we switch on the emission of neutrinos and antineutrinos at a constant rate
from a source located at $r_{0}$ moving in the positive outward direction. The
neutrinos are $99\%$ the electron flavor and $1\%$ $x$-flavor; the
antineutrinos are $90\%$ the electron antineutrino flavor and $10\%$ the
$x$-antineutrino flavor. For $r<c\,t$ the solution of the moment transport
equations is that

$\displaystyle\frac{E_{ee}^{\rm(th)}(r)}{{\rm Tr}[E^{\rm(th)}(r)]}$
$\displaystyle=[1-P^{(V)}_{T}(r)]\frac{E_{ee}(r_{0})}{{\rm
Tr}[E(r_{0})]}+P^{(V)}_{T}(r)\frac{E_{xx}(r_{0})}{{\rm Tr}[E(r_{0})]}$ (A1)
$\displaystyle\frac{E_{xx}^{\rm(th)}(r)}{{\rm Tr}[E^{\rm(th)}(r)]}$
$\displaystyle=[1-P^{(V)}_{T}(r)]\frac{E_{xx}(r_{0})}{{\rm
Tr}[E(r_{0})]}+P^{(V)}_{T}(r)\frac{E_{ee}(r_{0})}{{\rm Tr}[E(r_{0})]}$ (A2)

where the function $P^{(V)}_{T}$ is the flavor transition probability

$\displaystyle P^{(V)}_{T}(r)$
$\displaystyle=\sin^{2}(2\theta)\,\sin^{2}\left[\frac{\delta
m^{2}}{4p}(r-r_{0})\right],$ (A3)

$\theta$ is the mixing angle, $\delta m^{2}$ is the squared neutrino mass
difference, and $p$ is the neutrino energy. Analogous expressions to Eqs. (A1)
and (A2) for the theoretical anti-neutrino energy density moments exist with
an identical transition probability
$\overline{P}^{(V)}_{T}(r)=P^{(V)}_{T}(r)$.

Figure 9: Results from the vacuum test of FLASH plotted against radial
coordinate. [Top] Diagonal density matrix components of the energy density
moment ($E_{ii}$) normalized by the trace of $E$. Barred components denote the
anti-neutrino counterparts, i.e., $E_{\bar{e}\bar{e}}\equiv\overline{E}_{ee}$
and are scaled by the trace of $\overline{E}$. [Bottom] Errors in various
components. Green and black lines correspond to difference between the
calculated $E_{ii}$ and the theoretical predictions for Eqs. (A1) and (A2)
scaled by the trace of $E$. Anti-neutrino counterparts follow from similar
equations. Solid blue (dashed red) correspond to the difference between the
calculated length of the polarization vector and the initial length of the
polarization vector for neutrinos (anti-neutrinos) scaled by the initial
polarization vector length.

Figure 9 shows the results of this test at a time $t=3.3\times 10^{-5}\,{\rm
s}$. We use a mixing angle $\theta=0.28818$, a mass-squared difference $\delta
m^{2}=6.9\times 10^{-4}\,{\rm eV}^{2}$, and an energy $p=1\,{\rm MeV}$. The
top panel of the figure shows the flavor content evolving with the familiar
oscillatory pattern. We use the shorthand notation of
$E_{\bar{e}\bar{e}}\equiv\overline{E}_{ee}$ for the electron anti-neutrino
energy density moment, and similarly for the $\bar{x}\bar{x}$ component. In
the lower panel we plot the relative error of the numerical solution compared
to the analytic results in Eqs. (A1), (A2), and the anti-neutrino analogs,
namely

${\rm Error}=\frac{Q-Q^{{\rm(th)}}}{Q^{\rm(th)}},$ (A4)

for the relevant quantity $Q$. In addition to the errors in the energy density
moments, we give the relative error in the Bloch polarization vector for
neutrinos, $L$, and anti-neutrinos, $\overline{L}$. Notice that the largest
errors occur at small radii for $E_{ee}$ and $\overline{E}_{ee}$ (and
consequently also for $L$ and $\overline{L}$) at a level of a few parts in
$10^{3}$. This error is due to numerical oscillations at small radii which are
a result of the inner boundary condition in spherical coordinates and not
related to flavor-transformation. Overall, the FLASH calculation is able to
successfully reproduce the analytic results to a precision better than $\sim
1\%$.

###  A.2 MSW Oscillations

Our second 1D test case is similar to the vacuum test case described in the
previous section. In addition to the vacuum potential, we now add a constant
matter density to observe MSW oscillations (Bethe, 1986; Haxton, 1986). Our
domain ranges from the inner boundary at $r_{0}=10\,{\rm km}$ to the outer
boundary at $r=20\,{\rm km}$. We use 320 cells evenly-spaced in the radial
coordinate. The evolution equations for this test can actually be transformed
such that the analytic solution to the energy density moments is the same as
Eqs. (A1) and (A2) with a different MSW transition probability

$P^{(M)}_{T}(r)=\sin^{2}(2\theta_{\rm eff})\sin^{2}\left[\frac{\delta
m^{2}_{\rm eff}}{4p}(r-r_{0})\right],$ (A5) $\sin^{2}(2\theta_{\rm
eff})=\frac{\sin^{2}(2\theta)}{\sin^{2}(2\theta)+C^{2}},$ (A6) $\delta
m^{2}_{\rm eff}=\delta m^{2}\sqrt{\sin^{2}(2\theta)+C^{2}},$ (A7)
$C=\cos(2\theta)-\frac{2\sqrt{2}\,p\,G_{F}Y_{e}\rho}{\delta m^{2}\,m_{u}},$
(A8)

where $\rho$ is the matter density, $Y_{e}$ the electron fraction, and $m_{u}$
the atomic mass unit. For anti-neutrinos $\overline{P}^{(M)}_{T}$ is
calculated the same as above but with the placement of a minus sign on the
electron fraction, $Y_{e}\rightarrow-Y_{e}$. We can then compare this solution
to the FLASH calculation, which simply adds $H_{M}$ to $H_{V}$ in Eqs. (24)
and (25).

The results of this test case at $t=1.67\times 10^{-4}\,{\rm s}$ are shown in
Fig. 10, where the panels have the same plotting conventions as those in Fig.
9. We use a matter density of $\rho=8\times 10^{3}\,\mathrm{g\,cm}^{-3}$ and
an electron fraction $Y_{e}=0.5$. The familiar oscillatory evolution of the
neutrinos and anti-neutrinos as a function of distance is apparent. The
difference in period between the neutrinos and anti-neutrinos is due to the
replacement of $Y_{e}$ by $-Y_{e}$ in Eq. (A8). The errors for this MSW test
case are of a similar order of magnitude as the vacuum case.

Figure 10: Results from the MSW test of FLASH plotted against radial
coordinate. Notation for both panels is the same as Fig. 9. Error for the
$E_{ii}$ components corresponds to the difference between the calculated
$E_{ii}$ and the theoretical predictions using the transition probability in
Eq. (A5) for neutrinos, and a similar expression for anti-neutrinos.

###  A.3 Bipolar Oscillations

Our final test case is a study of bipolar neutrino oscillations. When using
negative energy eigenvalues for the anti-neutrino density matrices, systems
that are represented by two nearly oppositely directed flavor polarization
vectors are termed _bipolar systems_ (Duan et al., 2010). As they evolve, the
neutrino polarization vector nearly flips in direction, while the anti-
neutrino polarization vector exhibits the same behavior but $\pi$ radians out
of phase. In our treatment, we use positive energy eigenvalues for the anti-
neutrinos in Eq. (8), and so we begin with aligned flavor polarization vectors
to study the bipolar behavior. As a result, we would expect bipolar
oscillations in phase between the neutrinos and anti-neutrinos. We study a
simple case where the system is homogeneous and monochromatic. Only the
$H_{V}$ and $H_{\nu}$ terms are included in the Hamiltonian for Eqs. (24) and
(26). The flux moments are all zero. We pick $\delta m^{2}/p$ such that the
period of vacuum oscillations is $2\pi\,{\rm s}$. The energy density is picked
such that the strength of the self-interactions is $1\,{\rm Hz}$.

Figure 11 shows the results of our bipolar test for neutrinos — the anti-
neutrinos evolve in the same manner. The various color curves label the
individual components of the energy density matrix. Initially the neutrinos
are completely in $x$-flavor eigenstates. They subsequently flip to a state
which is $\sim 95\%$ electron flavor, before reverting back to the initial
state. This pattern repeats indefinitely as the system is always unstable.
Vertical dashed black lines give the locations in time of the period of
oscillations in $E_{xx}$ as predicted in Hannestad et al. (2006), centered on
the first peak of the $E_{ee}$ curve.333For the interested reader: Xiong et
al. (2023) gives an analytical formulation of the curves in Fig. 11. Note that
the period for the flavor off-diagonal components, ${\rm Re}(E_{ex})$ and
${\rm Im}(E_{ex})$ have double the period as the diagonal components. There
exists a slight decrease in the period of the FLASH results, which becomes
more apparent as time continues. At the first full period near $t\sim
4.1\,{\rm s}$, the relative difference in expected versus calculated time is
$5.7\times 10^{-3}$. At the second full period, near $t\sim 6.8\,{\rm s}$, the
drift is $6.9\times 10^{-3}$.

Figure 11: Results from the bipolar test of FLASH plotted against time. Only
shown are density matrix components for neutrinos. We include vertical dashed
black lines to indicate the period of the oscillations based on the analytical
result of Ref. Hannestad et al. (2006).

##  Appendix B Moments in an Orthornormal Tetrad

In order to perform simulations of realistic neutrino distributions in a small
domain using codes that assume a Minkowski metric, we need to transform the
radiation field quantities output by SpEC-Hydro into an orthonormal tetrad
comoving with the background fluid. The radiation field is given in terms of
the lab-frame energy density $E^{(\nu_{a})}$, energy flux $F_{i}^{(\nu_{a})}$,
the average neutrino energy $\epsilon^{(\nu_{a})}$, and the fluid transport
velocity $\widetilde{v}^{i}$ at each point in space and for each neutrino
flavor $a$. The spacetime metric is given in standard 3+1 quantities: the
lapse $\alpha$, shift $\beta^{i}$, and 3-metric $\gamma_{ij}$. This is a
common procedure within general-relativistic truncated moment simulation codes
(Shibata et al., 2011; Foucart et al., 2016b; Kuroda et al., 2016), but we
make the procedure explicit here for completeness. Note that in the main text,
we revert to the more conventional symbols $E$ and $\mathbf{F}$ to indicate
tetrad quantities, but they refer to the lab-frame quantities in this appendix
for consistency with Shibata et al. (2011).

The radiation stress-energy tensor can be defined in either the lab or
comoving frames as (Shibata et al., 2011)

$\displaystyle T^{\mu\nu}$
$\displaystyle=En^{\mu}n^{\nu}+F^{\mu}n^{\nu}+n^{\mu}F^{\nu}+P^{\mu\nu}\,\,\mathrm{[lab]}$
(B1)
$\displaystyle=Ju^{\mu}u^{\nu}+H^{\mu}u^{\nu}+u^{\mu}H^{\nu}+S^{\mu\nu}\,\,\mathrm{[comoving]}$

To get the moments in the comoving frame, we evaluate the full stress-energy
tensor and project onto suitably chosen tetrad basis vectors. This is made
more complicated by the fact that the pressure tensor needs to be expressed as
a linear combination of the optically thick and thin limits according to the
two-moment closure.

The index on the shift vector can be lowered with the three-metric
($\beta_{i}=\beta^{j}\gamma_{ij}$). The vector that is normal to the
hypersurface of constant time $t$ can be constructed as
$n_{\mu}=\alpha(-1,0,0,0)$. The transport velocity is defined in terms of the
four-velocity as $\widetilde{v}^{i}=u^{i}/u^{t}$, where $W=\alpha u^{t}$ is
the generalized Lorentz factor. From these it is straightforward to construct
the fluid four-velocity $u^{\alpha}$ and the three-velocity
$v^{i}=u^{i}/W+\beta^{i}/\alpha$.

The optically thin pressure tensor is straightforwardly defined as
$P^{ij}_{\mathrm{thin}}=(E/F^{k}F_{k})F^{i}F^{j}$, where contractions are
naturally evaluated using the lab-frame three-metric $\gamma_{ij}$. Obtaining
$P^{ij}_{\mathrm{thick}}$ is more involved, but can be determined analytically
as follows. It is straightforward to write the optically thick pressure tensor
in the lab frame as $S^{ij}_{\mathrm{thick}}=Jh^{ij}/3$, where
$h^{ij}=g^{ij}+u^{i}u^{j}$ is the comoving-frame three metric. By contracting
Eq. (B1) with various combinations of $u^{\alpha}$ and $n^{\alpha}$, one can
work out that the comoving-frame radiation energy density is [i.e., starting
with Eq. (87) in Foucart et al. (2016b)]

$\displaystyle J$ $\displaystyle=T^{\mu\nu}u_{\mu}u_{\nu}$ (B2)
$\displaystyle=W^{2}\left[E-2v^{i}F_{i}+JW^{2}(v_{i}v^{i})^{2}+2WH_{i}v^{i}(v_{i}v^{i})+S^{ij}v_{i}v_{j}\right]$

Plugging in $S^{ij}$ from above, this in turn yields the comoving-frame energy
density in the optically thick limit:

$J_{\mathrm{thick}}=\frac{3}{2W^{2}+1}\left[(2W^{2}-1)E-2W^{2}v^{i}F_{i}\right]$
(B3)

Similarly,

$\displaystyle H^{i}_{\mathrm{thick}}$
$\displaystyle=-T^{\mu\nu}u_{\mu}(g_{\nu}^{\phantom{\nu}i}+u_{\nu}u^{i})$ (B4)
$\displaystyle=\frac{F^{i}}{W}+\frac{v^{i}W}{2W^{2}+1}\left[(4W^{2}+1)F_{i}v^{i}-4W^{2}E\right]$

Finally, these can be used to calculate the pressure tensor in the optically-
thick limit as

$\displaystyle P^{ij}_{\mathrm{thick}}$
$\displaystyle=T^{\mu\nu}(g_{\mu}^{\phantom{\mu}i}+u_{\mu}u^{i})(g_{\mu}^{\phantom{\nu}j}+u_{\mu}u^{j})$
(B5)
$\displaystyle=\frac{J_{\mathrm{thick}}}{3}(4W^{2}v^{i}v^{j}+\gamma^{ij})+W(H^{i}_{\mathrm{thick}}v^{j}+v^{i}H^{j}_{\mathrm{thick}})$

By plugging $E$, $F^{i}$, and $P^{ij}$ [from Eq. (36)] into Eq. (B1), we now
have the complete stress-energy tensor $T^{\mu\nu}$.

For a comoving orthonormal tetrad, the timelike basis vector is the four
velocity ($\hat{t}_{\mu}=u_{\mu}$). We then choose a spacelike trial vector
$(0,1,0,0)$ (and similarly for the $y$ and $z$ trial vectors) and apply Gram-
Schmidt orthonormalization to get the spacelike basis vectors $x^{(i)}_{\mu}$.
The moments in the orthonormal tetrad are then

$\displaystyle J_{\mathrm{tet}}$
$\displaystyle=T^{\mu\nu}\hat{t}_{\mu}\hat{t}_{\nu},$ (B6) $\displaystyle
H^{i}_{\mathrm{tet}}$
$\displaystyle=T^{\mu\nu}\hat{t}_{\mu}\hat{x}^{(i)}_{\nu}.$ (B7)

From this point forward, we do not have to think about spacetime curvature, as
at each location spacetime is locally flat. Specifically, the moments can now
be used more intuitively. For example, the flux factor is simply
$f=\sqrt{\sum_{i}(H^{i}_{\mathrm{tet}})^{2}}/J_{\mathrm{tet}}$.

## References

  * Abbar (2023) Abbar, S. 2023, Phys. Rev. D, 107, 103006
  * Abbar & Duan (2018) Abbar, S., & Duan, H. 2018, Phys. Rev. D, 98, 043014
  * Abbar et al. (2019) Abbar, S., Duan, H., Sumiyoshi, K., Takiwaki, T., & Volpe, M. C. 2019, Phys. Rev. D, 100, 043004
  * Ardevol-Pulpillo et al. (2019) Ardevol-Pulpillo, R., Janka, H. T., Just, O., & Bauswein, A. 2019, Mon. Not. Roy. Astron. Soc., 485, 4754
  * Balantekin et al. (2023) Balantekin, A. B., Cervia, M. J., Patwardhan, A. V., Rrapaj, E., & Siwach, P. 2023, Eur. Phys. J. A, 59, 186
  * Balantekin & Pehlivan (2007) Balantekin, A. B., & Pehlivan, Y. 2007, J. Phys. G, 34, 47
  * Baxter et al. (2022) Baxter, A. L., et al. 2022, Astrophys. J., 925, 107
  * Bellm et al. (2018) Bellm, E. C., Kulkarni, S. R., Graham, M. J., et al. 2018, Publications of the Astronomical Society of the Pacific, 131, 018002
  * Bethe (1986) Bethe, H. A. 1986, Phys. Rev. Lett., 56, 1305
  * Bhattacharyya & Dasgupta (2020) Bhattacharyya, S., & Dasgupta, B. 2020, Phys. Rev. D, 102, 063018
  * Bhattacharyya & Dasgupta (2022) —. 2022, Phys. Rev. D, 106, 103039
  * Blaschke & Cirigliano (2016) Blaschke, D. N., & Cirigliano, V. 2016, Phys. Rev. D, 94, 033009
  * Bollig et al. (2021) Bollig, R., Yadav, N., Kresse, D., et al. 2021, ApJ, 915, 28
  * Bruenn et al. (2020) Bruenn, S. W., Blondin, J. M., Hix, W. R., et al. 2020, The Astrophysical Journal Supplement Series, 248, 11
  * Burrows et al. (2019) Burrows, A., Radice, D., & Vartanyan, D. 2019, MNRAS, 485, 3153
  * Cernohorsky & Bludman (1994) Cernohorsky, J., & Bludman, S. A. 1994, ApJ, 433, 250
  * Cervia et al. (2022) Cervia, M. J., Siwach, P., Patwardhan, A. V., et al. 2022, Phys. Rev. D, 105, 123025
  * Cirigliano et al. (2015) Cirigliano, V., Fuller, G. M., & Vlasenko, A. 2015, Physics Letters B, 747, 27
  * Cowperthwaite et al. (2017) Cowperthwaite, P. S., et al. 2017, Astrophys. J., 848, L17
  * Curtis et al. (2019) Curtis, S., Ebinger, K., Fröhlich, C., et al. 2019, ApJ, 870, 2
  * Curtis et al. (2023) Curtis, S., Miller, J. M., Fröhlich, C., et al. 2023, ApJ, 945, L13
  * Cusinato et al. (2022) Cusinato, M., Guercilena, F. M., Perego, A., et al. 2022, The European Physical Journal A, 58, 99
  * Dasgupta (2022) Dasgupta, B. 2022, Phys. Rev. Lett., 128, 081102
  * Dasgupta et al. (2009) Dasgupta, B., Dighe, A., Raffelt, G. G., & Smirnov, A. Y. 2009, Phys. Rev. Lett., 103, 051105
  * Dasgupta et al. (2017) Dasgupta, B., Mirizzi, A., & Sen, M. 2017, JCAP, 02, 019
  * Dasgupta et al. (2018) Dasgupta, B., Mirizzi, A., & Sen, M. 2018, Phys. Rev. D, 98, 103001
  * Duan et al. (2011) Duan, H., Friedland, A., McLaughlin, G., & Surman, R. 2011, J. Phys. G, 38, 035201
  * Duan et al. (2006) Duan, H., Fuller, G. M., & Qian, Y.-Z. 2006, Phys. Rev. D, 74, 123004
  * Duan et al. (2010) —. 2010, Annual Review of Nuclear and Particle Science, 60, 569
  * Dubey et al. (2009) Dubey, A., Antypas, K., Ganapathy, M. K., et al. 2009, Parallel Computing, 35, 512
  * Ehring et al. (2023a) Ehring, J., Abbar, S., Janka, H.-T., Raffelt, G., & Tamborra, I. 2023a, Phys. Rev. D, 107, 103034
  * Ehring et al. (2023b) —. 2023b, Phys. Rev. Lett., 131, 061401
  * Fernández et al. (2022) Fernández, R., Richers, S., Mulyk, N., & Fahlman, S. 2022, Phys. Rev. D, 106, 103003
  * Foucart et al. (2023) Foucart, F., Duez, M. D., Haas, R., et al. 2023, Phys. Rev. D, 107, 103055
  * Foucart et al. (2021) Foucart, F., Duez, M. D., Hebert, F., et al. 2021, Astrophys. J., 920, 82
  * Foucart et al. (2016a) Foucart, F., O’Connor, E., Roberts, L., et al. 2016a, Phys. Rev., D94, 123016
  * Foucart et al. (2016b) Foucart, F., Haas, R., Duez, M. D., et al. 2016b, Phys. Rev. D, 93, 044019
  * Freiburghaus et al. (1999) Freiburghaus, C., Rosswog, S., & Thielemann, F.-K. 1999, ApJ, 525, L121
  * Frensel et al. (2017) Frensel, M., Wu, M.-R., Volpe, C., & Perego, A. 2017, Phys. Rev. D, 95, 023011
  * Froustey et al. (2020) Froustey, J., Pitrou, C., & Volpe, M. C. 2020, JCAP, 12, 015
  * Froustey et al. (2023) Froustey, J., et al. 2023, In Preparation
  * Fryxell et al. (2000) Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273
  * Fujibayashi et al. (2023) Fujibayashi, S., Kiuchi, K., Wanajo, S., et al. 2023, Astrophys. J., 942, 39
  * George et al. (2022) George, M., Lin, C.-Y., Wu, M.-R., Liu, T. G., & Xiong, Z. 2022, arXiv:2203.12866, 2203.12866
  * Gizzi et al. (2019) Gizzi, D., O’Connor, E., Rosswog, S., et al. 2019, Mon. Not. Roy. Astron. Soc., 490, 4211
  * Grohs et al. (2016) Grohs, E., Fuller, G. M., Kishimoto, C. T., Paris, M. W., & Vlasenko, A. 2016, Phys. Rev. D, 93, 083522
  * Grohs et al. (2022) Grohs, E., Richers, S., Couch, S. M., et al. 2022, arXiv:2207.02214
  * Hannestad et al. (2006) Hannestad, S., Raffelt, G. G., Sigl, G., & Wong, Y. Y. Y. 2006, Phys. Rev. D, 74, 105010
  * Harten et al. (1983) Harten, A., Lax, P. D., & Leer, B. v. 1983, SIAM review, 25, 35
  * Haxton (1986) Haxton, W. C. 1986, Phys. Rev. Lett., 57, 1271
  * Hayashi et al. (2022) Hayashi, K., Fujibayashi, S., Kiuchi, K., et al. 2022, Phys. Rev. D, 106, 023008
  * Holmbeck et al. (2020) Holmbeck, E. M., et al. 2020, Astrophys. J. Suppl., 249, 30
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science Engineering, 9, 90
  * Izaguirre et al. (2017) Izaguirre, I., Raffelt, G., & Tamborra, I. 2017, Phys. Rev. Lett., 118, 021101
  * Johns (2023) Johns, L. 2023, Phys. Rev. Lett., 130, 191001
  * Johns & Nagakura (2021) Johns, L., & Nagakura, H. 2021, Phys. Rev. D, 103, 123012
  * Johns et al. (2020) Johns, L., Nagakura, H., Fuller, G. M., & Burrows, A. 2020, Phys. Rev. D, 102, 103017
  * Johns et al. (2020) Johns, L., Nagakura, H., Fuller, G. M., & Burrows, A. 2020, Phys. Rev. D, 101, 043009
  * Johns & Xiong (2022) Johns, L., & Xiong, Z. 2022, Phys. Rev. D, 106, 103029
  * Jones et al. (2001) Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python, , , [Online; accessed <today>]
  * Just et al. (2022) Just, O., Abbar, S., Wu, M.-R., et al. 2022, Phys. Rev. D, 105, 083024
  * Just et al. (2015) Just, O., Obergaulinger, M., & Janka, H.-T. 2015, Monthly Notices of the Royal Astronomical Society, 453, 3386
  * Kalogera et al. (2021) Kalogera, V., et al. 2021, arXiv:2111.06990
  * Kneller et al. (2023) Kneller, J., et al. 2023, In Preparation
  * Kuroda et al. (2016) Kuroda, T., Takiwaki, T., & Kotake, K. 2016, The Astrophysical Journal Supplement Series, 222, 20
  * Lacroix et al. (2022) Lacroix, D., Balantekin, A. B., Cervia, M. J., Patwardhan, A. V., & Siwach, P. 2022, Phys. Rev. D, 106, 123006
  * Lien & Fields (2009) Lien, A., & Fields, B. D. 2009, J. Cosmology Astropart. Phys, 2009, 047
  * Lippuner & Roberts (2015) Lippuner, J., & Roberts, L. F. 2015, ApJ, 815, 82
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration, Abell, P. A., Allison, J., et al. 2009, arXiv e-prints, arXiv:0912.0201
  * Malkus et al. (2014) Malkus, A., Friedland, A., & McLaughlin, G. C. 2014, arXiv:1403.5797
  * Malkus et al. (2012) Malkus, A., Kneller, J. P., McLaughlin, G. C., & Surman, R. 2012, Phys. Rev. D, 86, 085015
  * Martin et al. (2023a) Martin, J. D., Neill, D., Roggero, A., Duan, H., & Carlson, J. 2023a, arXiv e-prints, arXiv:2307.16793
  * Martin et al. (2023b) Martin, J. D., Roggero, A., Duan, H., & Carlson, J. 2023b, arXiv e-prints, arXiv:2301.07049
  * Martin et al. (2020) Martin, J. D., Yi, C., & Duan, H. 2020, Physics Letters B, 800, 135088
  * McLaughlin et al. (1996) McLaughlin, G. C., Fuller, G. M., & Wilson, J. R. 1996, Astrophys. J., 472, 440
  * Meurer et al. (2017) Meurer, A., Smith, C. P., Paprocki, M., et al. 2017, PeerJ Computer Science, 3, e103
  * Miller et al. (2020) Miller, J. M., Sprouse, T. M., Fryer, C. L., et al. 2020, Astrophys. J., 902, 66
  * Minerbo (1978) Minerbo, G. N. 1978, J. Quant. Spec. Radiat. Transf., 20, 541
  * Morinaga (2022) Morinaga, T. 2022, Phys. Rev. D, 105, L101301
  * Mukhopadhyay (2022) Mukhopadhyay, P. 2022, PhD thesis, Stanford U.
  * Myers et al. (2022) Myers, M., Cooper, T., Warren, M., et al. 2022, Phys. Rev. D, 105, 123036
  * Nagakura (2022) Nagakura, H. 2022, Phys. Rev. D, 106, 063011
  * Nagakura (2023) Nagakura, H. 2023, Phys. Rev. Lett., 130, 211401
  * Nagakura et al. (2021) Nagakura, H., Johns, L., Burrows, A., & Fuller, G. M. 2021, Phys. Rev. D, 104, 083025
  * Neilsen et al. (2014) Neilsen, D., Liebling, S. L., Anderson, M., et al. 2014, Phys. Rev. D, 89, 104029
  * O’Connor & Couch (2018) O’Connor, E. P., & Couch, S. M. 2018, ApJ, 854, 63
  * Padilla-Gay et al. (2021) Padilla-Gay, I., Shalgar, S., & Tamborra, I. 2021, JCAP, 01, 017
  * Patwardhan et al. (2021) Patwardhan, A. V., Cervia, M. J., & Balantekin, A. B. 2021, Phys. Rev. D, 104, 123035
  * Pehlivan et al. (2011) Pehlivan, Y., Balantekin, A. B., Kajino, T., & Yoshida, T. 2011, Phys. Rev. D, 84, 065008
  * Perego et al. (2016) Perego, A., Cabezón, R., & Käppeli, R. 2016, Astrophys. J. Suppl., 223, 22
  * Press et al. (1993) Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1993, Numerical Recipes in FORTRAN; The Art of Scientific Computing, 2nd edn. (New York, NY, USA: Cambridge University Press)
  * Radice et al. (2022) Radice, D., Bernuzzi, S., Perego, A., & Haas, R. 2022, Mon. Not. Roy. Astron. Soc., 512, 1499
  * Reichert et al. (2021) Reichert, M., Obergaulinger, M., Eichler, M., Aloy, M. Á., & Arcones, A. 2021, MNRAS, 501, 5733
  * Richers (2020) Richers, S. 2020, Phys. Rev. D, 102, 083017
  * Richers (2022) —. 2022, Phys. Rev. D, 106, 083005
  * Richers et al. (2021a) Richers, S., Willcox, D., & Ford, N. 2021a, Phys. Rev. D, 104, 103023
  * Richers et al. (2021b) Richers, S., Willcox, D. E., Ford, N. M., & Myers, A. 2021b, Phys. Rev. D, 103, 083013
  * Richers et al. (2019) Richers, S. A., McLaughlin, G. C., Kneller, J. P., & Vlasenko, A. 2019, Phys. Rev. D, 99, 123014
  * Roggero et al. (2022) Roggero, A., Rrapaj, E., & Xiong, Z. 2022, Phys. Rev. D, 106, 043022
  * Rosswog & Liebendörfer (2003) Rosswog, S., & Liebendörfer, M. 2003, MNRAS, 342, 673
  * Rrapaj (2020) Rrapaj, E. 2020, Phys. Rev. C, 101, 065805
  * Ruffert & Janka (1999) Ruffert, M., & Janka, H.-T. 1999, Astron. Astroph., 344, 573
  * Sawyer (2005) Sawyer, R. F. 2005, Phys. Rev. D, 72, 045003
  * Shibata et al. (2011) Shibata, M., Kiuchi, K., Sekiguchi, Y.-i., & Suwa, Y. 2011, Progress of Theoretical Physics, 125, 1255
  * Sigl & Raffelt (1993) Sigl, G., & Raffelt, G. 1993, Nuclear Physics B, 406, 423
  * Siwach et al. (2023) Siwach, P., Suliga, A. M., & Balantekin, A. B. 2023, Phys. Rev. D, 107, 023019
  * Skinner et al. (2019) Skinner, M. A., Dolence, J. C., Burrows, A., Radice, D., & Vartanyan, D. 2019, The Astrophysical Journal Supplement Series, 241, 7
  * Stapleford et al. (2020) Stapleford, C. J., Fröhlich, C., & Kneller, J. P. 2020, Phys. Rev. D, 102, 081301
  * Surman et al. (2006) Surman, R., McLaughlin, G. C., & Hix, W. R. 2006, Astrophys. J., 643, 1057
  * Tartaglia et al. (2018) Tartaglia, L., et al. 2018, Astrophys. J., 853, 62
  * Tian et al. (2017) Tian, J. Y., Patwardhan, A. V., & Fuller, G. M. 2017, Phys. Rev. D, 96, 043001
  * Tian et al. (2017) Tian, J. Y., Patwardhan, A. V., & Fuller, G. M. 2017, Phys. Rev. D, 95, 063004
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science Engineering, 13, 22
  * Vlasenko et al. (2014) Vlasenko, A., Fuller, G. M., & Cirigliano, V. 2014, Phys. Rev. D, 89, 105004
  * Vlasenko & McLaughlin (2018) Vlasenko, A., & McLaughlin, G. C. 2018, Phys. Rev. D, 97, 083011
  * Volpe et al. (2013) Volpe, C., Väänänen, D., & Espinoza, C. 2013, Phys. Rev. D, 87, 113010
  * Volpe (2023) Volpe, M. C. 2023, arXiv:2301.11814
  * Wanajo et al. (2014) Wanajo, S., Sekiguchi, Y., Nishimura, N., et al. 2014, Astrophys.J.Lett., 789, L39
  * Willcox & Richers (2021) Willcox, D. E., & Richers, S. 2021, AMReX-Astro/Emu:, v.v1.3, Zenodo, doi:10.5281/zenodo.5514018
  * Wu et al. (2016) Wu, M.-R., Duan, H., & Qian, Y.-Z. 2016, Phys. Lett. B, 752, 89
  * Wu & Tamborra (2017) Wu, M.-R., & Tamborra, I. 2017, Phys. Rev. D, 95, 103007
  * Xiong et al. (2022) Xiong, Z., Johns, L., Wu, M.-R., & Duan, H. 2022, arXiv:2212.03750
  * Xiong et al. (2023) Xiong, Z., Wu, M.-R., Martínez-Pinedo, G., et al. 2023, Phys. Rev. D, 107, 083016
  * Xiong et al. (2023) Xiong, Z., Wu, M.-R., & Qian, Y.-Z. 2023, Phys. Rev. D, 108, 043007
  * Zaizen & Morinaga (2021) Zaizen, M., & Morinaga, T. 2021, Phys. Rev. D, 104, 083035
  * Zhang & Burrows (2013) Zhang, Y., & Burrows, A. 2013, Phys. Rev. D, 88, 105009