Two-moment Neutrino Flavor Transformation with Applications to the Fast Flavor Instability in Neutron Star Mergers

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 the flavor-transformation with spatial transport to evolve the neutrino radiation quantum field. We show that moments capture the dynamical features of fast flavor instabilities in a variety of systems, although our technique is by no means a universal blueprint for solving fast flavor transformation. To evaluate the effectiveness of our moment results, we compare to a more precise particle-in-cell method. Based on our results, we propose areas for improvement and application to complementary techniques in the future.


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 neu-tron 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 CCSNe simulations (Stapleford et al. 2020;Ehring et al. 2023a,b); and the results of element synthesis obtained in neutrino-cooled accretion disks (e.g., Malkus et al. 2012;Just et al. 2022), supernovae (e.g., Duan et al. 2011;Mukhopadhyay 2022;Fujimoto & Nagakura 2023) and hypermassive neutron star outflows (e.g., Fernández et al. 2022;George et al. 2020;Li & Siegel 2021).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 neu-trino 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 (e.g., Nagakura et al. (2014Nagakura et al. ( , 2018)); Iwakami et al. (2020) for CCSNe and Miller et al. (2019); Miller et al. (2019) for NSMs).
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(Malkus et al. , 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 (Dasgupta et al. 2018;Johns & Nagakura 2021;Nagakura & Johns 2021;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 instabili-ties (Johns 2023;Johns & Xiong 2022;Xiong et al. 2023;Xiong et al. 2023a).
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;Johns et al. 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. 2023).
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 ℏ = c = 1.When giving results of numerical calculations, we use cgs units.

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, ϱ = ϱ(t, x, p), and corresponding generalized density matrix ϱ for antineutrinos, which are functions of time t, spatial location x, and momentum 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 ϱ and ϱ simply "density matrices" for the sake of brevity.If we are describing neutrinos and anti-neutrinos with 2 chiral states, ϱ and ϱ are 2 n f × 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 × n f for each of the neutrinos and anti-neutrinos.In the case of three flavors of neutrinos, namely e, µ, τ , we can write the neutrino density matrix as the following 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 appropri-ate phase-space integrals of ϱ: where the superscript i indicates a component of a 3vector.Note that we approximate neutrinos as ultrarelativistic by setting the neutrino energy equal to the 3-momentum magnitude p and that N , E, and the components of F are all n f × n f matrices.The expressions for anti-neutrinos are analogous.
The QKEs for ϱ and ϱ are where the single dot over a variable indicates differentiation with respect to time.In Eqs. ( 5) and ( 6), C and 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, denoting the vacuum (H V ), matter (H M ), and selfinteraction (H ν ) terms, and where * denotes complex conjugation.The vacuum term arises from non-zero neutrino rest masses, and we write it as where U is the PMNS matrix and 3 ) 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 in the case the matter fluid has zero velocity where G F ≃ 1.166 × 10 −11 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 = diag(1, 0, 0) for three flavors.We will work in a frame comoving with the matter fluid implying Eq. ( 10) is valid (see App. B).
Finally, the self-interaction potential is a consequence of neutrinos interacting with the background of other neutrinos (11) where ϑ is the angle between the free variable p and the integration variable q.

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 ϱ to be where i, j ∈ {x, y, z} are spatial indices.Note that a different convention was chosen with respect to Myers et al. (2022) (no 1/4π prefactors), for consistency with Grohs et al. (2023).Analogous expressions exist for the antineutrinos.The integrals in these definitions are only over the momentum-space solid angle Ω 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, x, p) via 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.
In this manuscript, we will only consider mono-energetic neutrinos, and as a result, our expressions for number density and number density moment differ by a factor of energy-bin width ∆p.Note that F is the specific energy flux, but 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 with the norm for a component of the flavor matrix defined as where i runs over the spatial indices x, y, z.
From comparing the definitions of the moments and Eq. ( 11) for self-interactions, we observe that the selfinteraction term can be written as where the moment-self-interaction terms are We will use these moment self-interaction expressions for evolving our dependent variables of E and F j .Note, however, that when writing Eqs. ( 22) and ( 23), the energy-integrated quantities N and F j appear.Indeed, for the particular physical phenomena we study herenamely the FFI -the ELN crossing depends on the number moments and not the energy ones (E and intensity).The simulations we present in this work are for mono-energetic neutrino distributions where N and E are equal up to a units factor as shown in Eq. ( 15).Nevertheless, we make the distinction between N and E for the FFI under the guise of an eventual incorporation of multi-energy distributions.
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 µν = 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.In addition, we will assume that gradients of the fluid velocity are locally approximately zero, which further simplifies the equations of motion.
With everything defined we can now write out the evolution equations for E and F by performing moment integrations of Eqs. ( 5) and ( 6) and scaling by appropriate factors of p 3 /(2π) 3 .In Cartesian coordinates the transport equations for the neutrino moments are where we have ignored the collision and force terms and assume the background matter to be homogeneous with zero velocity in the tetrad frame.Note that this form would also be applicable if we assume a Minkowski metric with nonzero velocity gradients for energy-integrated moments, but here we explicitly assume zero velocity everywhere.Note also 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 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 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 in the laboratory (Euler) frame where J α (e) = u α (n e − ne )I e .In the above Eq.( 29), F T denotes the transpose of the row vector F into a column vector.In addition, we define the four-velocity of the reference frame u µ = (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 In Eqs. ( 32) and (33), we adopt the convention where repeated indices are contracted with respect to the metric, i.e., A α A α = g αβ A α A β with a (−, +, +, +) convention.

The Maximum Entropy Closure
The evolution equations for the fluxes F and F involve the spatial gradients of the pressure tensors P and 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 momentumspace such that the angular entropy is extremized (Minerbo 1978;Cernohorsky & Bludman 1994).In other words, the neutrinos are distributed in the momentumspace 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, ψ aa , is (Minerbo 1978;Cernohorsky & Bludman 1994) where E aa is the energy density moment for species a, and µ = Ftet • Ω gives the angular dependence.Ω is the direction unit vector and F tet = F aa /F aa .The parameter Z aa follows from the constraint on the magnitude of the flux factor vector f 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 twomoment 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 where the thin and thick limits are 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 χ: 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.

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 antineutrino distributions) herald flavor instabilities (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.We note that although only an ELN crossing is initially present, an XLN crossing can subsequently appear during the evolution.
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.
We will determine whether an ELN crossing exists using energy density distributions, yet an ELN crossing utilizes number density distributions by definition.However, we stress that our energy and number variables are simply related by a constant of proportionality for mono-energetic distributions, and as a result, we will continue using the energy quantities below.Let E ee and Z ee define the maximum entropy distributions for electron neutrinos, and similarly E ee and Z ee for electron anti-neutrinos [see Eqs. ( 34) and ( 35)].Furthermore, we assume that the flux factors are separated by an angle θ, i.e., The distributions cross if (Richers 2022) where 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.

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.

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: ν e , ν e , and ν other for all other neutrinos.We modify the classical code by by distinguishing ν other from ν 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 µµ , Re(E eµ ), and Im(E eµ ) and the four charge-conjugate counterparts for the antineutrinos.In reality, neutrinos would oscillate between e and the other two flavors, namely, µ and τ .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 cells and group cells together into blocks to parallelize the computation over processors.Each block contains 16 3 cells along with ghost cells.The choice of 16 3 cells per block results, in part, from the computational resources we use for this work.For a different platform, we would be free to change the size of the blocks depending on the number of cores and available memory.We use 6 ghost cells in each dimension so communication occurs only at the end of each full timestep, given a stencil size of 2 in each direction and 3 substeps within each full step1 .The threestep 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 timestep 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 × 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 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 timestep 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 timestep 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 offdiagonal components of E are zero.We would also need to apply the same unitary transformation to each vector component of 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 χ for neutrinos and a single χ for antineutrinos using flavor-traced flux factors.Specifically, those flavor-traced flux factors are defined as and a similar expression for the anti-neutrinos and f (F T ) , where xi are the Cartesian unit vectors.These flavor traced flux factors are substituted into Eq.( 39) to obtain χ and χ, 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., ϱ ee , ϱ xx , Re[ϱ ex ], and Im[ϱ ex ].In other words, Eqs. ( 36), (37), and (38) all have flavor indices on each quantity, except for χ.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 ∼ f xx ∼ f (F T ) and similarly for the antineutrino flux factors.We leave a more detailed analysis of possible closures to future work (Kneller et al. 2023).

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 fourthorder 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).

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 ϱ xx and ϱ xx , since all heavy lepton neutrino and antineutrino species are gathered in a single species ν 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 selfinteracting 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 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 the 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 where −1 < A, B < 1 are uniform random numbers at each location x, N cc are the initial number density moments for the diagonal components, and the Hermitian perturbation is only applied for a ̸ = 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 implying the initial perturbations for the off-diagonal components of E and F are correlated.Analogous expressions exists for the anti-neutrino moments.
In the EMU calculations, each particle is assigned a 4momentum vector, weight, and density matrix.The 4momentum 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 δE ab to be reduced by √ 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.We give more details on the differences in the initial conditions between FLASH and EMU in App. C. 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.

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 Im(Ω) to denote the growth rate of |N ex | during instability.In addition, we use the symbol |k| max to denote the fastest growing mode in the discrete Fourier transform of N ex during instability.

3D Test Problems
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 solutions2 , 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 flavordiagonal number density moment.All three tests start with non-zero numbers of electron neutrinos and antineutrinos, 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).The MEC distributions are 3D as emphasized above and as a result, we plot polar representations of 2D crosssectional slices in the top row of Fig. 1.We measure the polar angle ϑ ∈ [0, 2π] counter-clockwise from the z axis.We take the slices such that the maximum values of the distributions (that is the directions of the fluxes F ee and F ee ) 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

Name
Nee N ee ΣN (x) fee f ee f (x) L Ngp (10 32 cm −3 ) (10 32 cm −3 ) (10 32 cm −3 ) (cm) Fiducial 4.89 4.89 0 (0, 0 , 1/3 ) (0, 0 , −1/3 ) (0, 0, 0)  apparent in the standard plots on the lower panel.The vectors originating from the origin on the polar plots show the peak direction of the distributions.The difference of the blue and red vectors is shown in dashed purple -for instance, it is coincident with the blue vector on the Fiducial case, and the vector difference is shown vividly in the 90Degree test.For all three tests, the coordinates are chosen so that the lepton number flux (i.e., the purple vector) lies along the z axis.We then orient the polar plane so that this axis points in the rightwards direction, and indicate this direction as ϑ = 0 in the lower plots.As in the polar plots, the blue and red curves give the electron neutrino and the electron anti-neutrino distributions.Here, the purple curve 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.The horizontal axis is t − tsat where tsat differs between test case and method of calculation.We use this definition for visualization purposes and stress the calculations are not simultaneous (see Fig. 15 for the same plots using directly the simulation time for the horizontal axis).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 flavormatrix N (i.e., number density of νe) scaled by Nee(t = 0).[Bottom] Plotted is the magnitude of the off-diagonal component of N scaled by Nee(0).For 3-flavor EMU calculations, we take the eµ component of N .
flavor EMU, and 3-flavor EMU simulations.We emphasize that, given the different initial perturbations for FLASH and EMU calculations (see Section 4.3), we do not expect identical time-evolution.Subsequently, the saturation time, t sat , when the off-diagonal terms saturate (located at the peak of the |N ex |/N ee (0) curves) depends on the initial conditions, as well as the kind of calculation.To aid in visualization when comparing the growth, saturation, and decoherence phases between FLASH and EMU, we define the horizontal axes in Fig. 2 as t − t sat using a different t sat for each calculation.We stress that none of the calculations are simultaneous with one another in simulation time -the alignment at t − t sat = 0 is a construct of the plot.Finally, 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 ⟨|N ex |⟩, evidencing one of the defining characteristics of the FFI.⟨|N ex |⟩ 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 | ≲ 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 ⟨N ee ⟩ 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 Im(Ω).Specifically for the Fiducial test case, Im(Ω) = 7.1 × 10 10 s −1 in FLASH, as compared to Im(Ω) = 6.3 × 10 10 s −1 in EMU, for a difference of ∼ 10%.Results are similar for the 90Degree test, with Im(Ω) = 5.4 (4.4) × 10 10 s −1 in FLASH (EMU).
Also, the asymptotic values for ⟨N ee ⟩ are nearly the same, with differences of ∼ 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.Note that during the post-saturation decoherence, there appear to be two phases indicated by different slopes in N ex .In the first, immediately after saturation, N ex decreases rapidly, but then numerical artifacts take over and decrease the decoherence rate (e.g., at around t − t sat = 0.75 ns in the Fiducial case).In the case of the FLASH calculations, this is due to the numerical diffusion from finite grid spacing, and in the case of EMU this is due to the finite number of computational particles that achieve a state of random uncorrelated fluctuations, the amplitude of which scale very slowly as N −1/2 p .Therefore the decoherence phase right after saturation is a robust physical prediction, but the late-time values of N ex show numerical artifacts.
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 ⟨N ee ⟩.The growth rate is faster for FLASH by ∼ 40% and the loss of coherence falls off faster.There is a smaller amount of time when |N ex | ≲ N ee , and thus less oscillations in the flavor-diagonal term.The result is an asymptotic value which is ≳ 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 ⟨N ee ⟩ 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 ⟨N ee ⟩ 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 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 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 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, canceling the omnipresent ELN one (Nagakura & Zaizen 2022;Zaizen & Nagakura 2023;Xiong et al. 2023b).Furthermore, the ν x and ν 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 ⟨N ee ⟩ 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 nor 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).
Pressure Moment -As discussed above, the MEC need not be a true representation of the distribution even if we find flavor-convergence in the asymptotic limit.Figure 3 gives a plot of the zz pressure tensor component for the electron neutrinos in the Fiducial test case.We pick the zz component for P ee as z is the direction of the net neutrino flux.For the geometry of the Fiducial test case, the thin and thick components of the pressure tensor reduce simply to P zz thin /E = 1 and P zz thick /E = 1/3, implying that the interpolated value from Eq. ( 36) is P zz /E = χ.For the purposes of analyzing our moment and PIC simulations, we plot the averaged values of P zz ee /E ee ∼ χ against the time as measured from the saturation peak.We choose this representation of our data as we do not expect qualitative differences for different cells.In contrast, Nagakura & Zaizen (2023) plots the time-averaged values of the pressure against the radial coordinate when comparing multi-angle results to closure approximations in a global CCSN simulation, showing the transition from the optically thick to thin limit.The solid black curve in Fig. 3 corresponds to the baseline EMU calculation, i.e., the solid black curve in the upper-left plot of Fig. 2. To calculate P zz ee for EMU, we use Eq. ( 14) to sum over the particles in a given cell and obtain the second angular moment of the distribution.We subsequently average over the simulation domain and normalize by the energy density moment.The dashed red curve gives the same quantity for the FLASH simulation.For FLASH, we first calculate the domain average of the energy and flux density moments.Along with the flavor-traced Eddington factor from Eq. ( 46), we use Eq. ( 36) with the averaged E ee and F ee to obtain P zz ee .Finally, we normalize by ⟨E ee ⟩.The constant value of the red dashed line shows that FLASH is conserving both neutrino energy density (i.e., particle number) and neutrino flux density (i.e., 3-momentum).
For diagnostic purposes, we include two other pressure quantities in Fig. 3.The dashed orange curve gives the pressure using the output FLASH energy and flux moments along with the classical MEC prescription (i.e., an Eddington factor calculated using pure diagonal fluxfactors without a trace over flavor).The solid blue curve gives the same but for output EMU quantities.By comparing the blue and black curves, we see how much the distribution in the PIC calculation differs from the classical MEC.Note that at times t < t sat , a finite number of particles causes the black curve to deviate from the blue one (we verified that increasing the number of particles reduces this discrepancy).After saturation, the black curve exhibits a larger amplitude of oscillations as compared to the blue curve.For large values of ⟨P zz ee ⟩, the actual PIC calculation is more forward peaked than the MEC approximation.The opposite would be true for small values of ⟨P zz ee ⟩, although it appears that the black and blue curves do not differ much at their minima.This finding is consistent with Nagakura & Zaizen (2023) during periods of significant flavor transformation [see Figs. ( 8) and ( 9) of the aforementioned work], despite the differences in plotting axes.Nagakura & Za-izen (2023) show the closure relation cannot always describe the shape of the flavor-transformed distribution during rapid flavor oscillations as the Eddington factor falls outside of the classically allowed range.Although ⟨P zz ee ⟩/⟨E ee ⟩ always falls within the classically allowed range for our Fiducial case, the difference between the blue and black lines is most acute at the maxima, and implies the MEC does not capture the multi-angle distribution at all times.Lastly, the black and blue curves oscillate nearly in phase with one another, indicating that the MEC contains the correct scaling of P zz ee with E ee and F ee but not the correct sensitivity.
We notice another difference in sensitivity when comparing the solid blue curve to the dashed orange curve of FLASH.The classical MEC calculation using FLASH data shows a smaller amplitude of oscillation, along with a larger frequency.We attribute the smaller amplitudes to the fact that the MEC underestimates the degree of forward-peaking of the distribution.The larger frequency correlates with the smaller time scales exhibited by FLASH, and observed in all three test cases.Notice that the dashed orange and solid blue curves do asymptote to similar values during the decoherence period, implying that the zeroth and first moments have similar values between the two methods of calculation.Finally, we note that our choice of utilizing the flavor-traced Eddington factor (dashed red curve in Fig. 3) over the classical MEC in the FLASH simulation results in a value of ⟨P zz ee ⟩ differing by ∼ 1% of ⟨E ee ⟩.As we operate in the optically thick limit at all times for this test case, we do not foresee that adopting the classical MEC prescription for calculating χ would alter the results in Fig. 2 by more than a few percent.
Fourier space analysis -We have discussed averages of the number density and pressure moments when presenting Figs. 2 and 3.In Fig. 4 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 0.00 0.25 0.50 0.75 1.00 1.25 1.50 t − t sat (10 −9 s) gives the same for the FLASH simulation.Also included are diagnostic quantities for FLASH (dashed orange) and EMU (solid blue) using a classical MEC along with the number and flux moments as given by the simulations.All quantities are averaged over the simulation domain.
of grid points per cm and a smaller maximum value of k.
The three DFTs in each panel of Fig. 4 are all from a time before saturation during the growth period: ∼ 0.1 ns before saturation for the Fiducial and 90Degree cases; ∼ 0.35 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 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. 4 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. 4. Soon after the simulations begin the DFT exhibits a peak with an associated wave number as evidenced in Fig. 4. 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| 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| max = 3.9 (3.1) cm −1 for FLASH (EMU) in the Fiducial test case, and |k| max = 3.1 (2.4) 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.
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 four are deduced from Fig. 2 and are the following, respectively: the maximum value of ⟨|N ex |/N ee (0)⟩ in the bottom panels; the ratio of ⟨|N ex |⟩ at the saturation time and a time t dec = t sat + 0.2 ns during the decoherence phase; the asymptote of ⟨N ee /N ee (0)⟩ in the top panels; and the slope of the line (in semi-log space) in the bottom panels.The fifth column gives the value of k at the peak of the DFT in Fig. 4. We give an uncertainty in parentheses for |k| max due to the finite box size, namely, δk = ±π/L.All tests show the FLASH calculations have larger values of Im(Ω) and |k| max compared to EMU.In addition, the rate of decline of |N ex | is larger for the FLASH calculations in all three tests.However, even with the different growth and loss of coherence rates, oscillations occur while the average value of |N ex | exhibits 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.Moreover, it varies with the random initial perturbations and should not be taken as a robust prediction for each calculation, contrary to the growth rate, instability lengthscale, and amount of flavor transformation.

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.Table 2. Numerical results for FLASH and EMU (2f) calculations for the three 3D test simulations.First column gives the ratio ⟨|Nex|/Nee(0)⟩ when t = tsat.Second column gives the ratio of off-diagonal magnitudes at two different times: tsat and t dec = tsat + 0.2 ns.Third column gives the asymptotic ratio ⟨Nee/Nee(0)⟩ post-saturation.Fourth column gives the growth rate Im(Ω) when the system is unstable in units of 10 10 s −1 .Last column gives the fastest growing mode in the domain |k|max in units of cm −1 , with an associated uncertainty in parentheses.The values of the first two columns are much more variable with the initial random perturbations than the last three.
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.
We analyze and simulate neutrino distributions at a selection of points in the polar slice of a snapshot at 5 ms post merger shown in Fig. 5. Green contours give matter densities of {10 11 , 10 12 , 10 13 , 10 14 } 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. (2023).
Table 3 and Fig. 6 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. 6 roughly correspond to the ẑ direction in Fig. 5.The leftward direction in the polar plot for point 3 in Fig. 6 roughly corresponds to the x direction in Fig. 5.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. 5.
Time evolution and FFI -Figure 7 shows the timeevolution 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. 7.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µ component in the bottom panel.
In all three NSM points, we see growth, saturation, and decoherence phases as we did Sec.5.1 and Fig. 2. Growth of |N ex | begins soon after the start of the sim-ulation and proceeds until ⟨|N ex |⟩ ∼ 0.1⟨N ee ⟩ in both FLASH and EMU calculations.Rapid oscillations develop and effect a decrease in ⟨N ee ⟩ towards complete flavormixing.Notice that for the first NSM point treated by FLASH, the conditions are such that ⟨N ee ⟩ 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) ⟨N ee ⟩ 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 ⟨N ee ⟩ 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), Im(Ω) = 5.2 (1.1) × 10 10 s −1 for this point.The bottom panels of Fig. 6 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.
We show in Fig. 8 3D-volume-renderings of the FLASH simulation for the first NSM point at four different simulation times: t = 0.01 ns; t = 0.18 ns; t = 0.21 ns; and t = 0.59 ns.The contours in each panel are for the phase, (10 32 cm −3 ) (10 32 cm −3 ) (10 32 cm −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 Ngp = 128 3 grid cells due to the longer wavelength of the fastest growing mode.ϕ ex , of the complex number N ex .The spatial structure in ϕ 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: ϕ ex = −π/2 (blue); ϕ ex = 0 (white); and ϕ ex = π/2 (red).The first panel shows a time close to the start of the simulation, where little flavortransformation 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 wave-0.60.8 1.0  16 for the same plots using directly the simulation time for the horizontal axis).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µ.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.
length 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 ϕ ex .The pattern seen here in Fig. 8 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 ϕ ex again shows that the moment method reproduces many features of the FFI on large and small scales.
Pressure Moment -In analogy to the Fiducial test case in Sec.5.1, we compare components of the pressure tensor between FLASH and EMU calculations for the NSM 1 point in Fig. 9.For the Fiducial test case, we chose the pressure component along the symmetry axis z when drawing Fig. 3.No such symmetry axis exists for the NSM 1 point, so we instead rotate into a primed reference frame where a principal axis is aligned with the flux vector for a given density matrix component.In other words, if we define a basis ( x ′ , y ′ , z ′ ) such that z ′ is the unit vector in the F aa direction, we compute P [rot] aa ≡ P z ′ z ′ aa , which we obtain after an appropriate spatial rotation of the pressure tensor.The rotation is different for each flavor, i.e., a = e, x, and computed at each time step.The solid black curve represents this quantity, averaged over the simulation domain and normalized by ⟨E aa ⟩, for electron and heavy lepton flavor neutrinos.If the pressure moment is obtained from the closure relation (36), then by construction P

[rot]
aa /E aa = χ aa .This is indeed the case for the dashed red curve in the top and bottom panels of Fig. 9, which is equal to the flavor-traced Eddington factor used in FLASH.For diagnostic purposes, we also represent the pressure moment computed using the classical MEC prescription, i.e., the non flavor-traced Eddington factor obtained from the first two angular moments in EMU (solid blue curve) and FLASH (dashed orange curve).
Similar to the Fiducial test case and Fig. 3, the MEC is able to capture some of the features of the underlying distribution.We note that the solid blue curve tracks the black curve quite closely and asymptotes to nearly identical values for the electron flavor pressure (similar to Fig. 3, the initial discrepancy between the black and blue curves is due to the finite number of particles).Similar to the Fiducial test case, the black curve tends to have more extreme maxima and minima, implying that the MEC underestimates the degree of forward-peaking (Nagakura & Zaizen 2023).For FLASH, the dashed orange curve follows the solid blue curve for roughly half a period during the onset of rapid flavor oscillations immediately after saturation.These oscillations terminate prematurely for the FLASH simulation and continue for EMU, implying final asymptotic values for ⟨P [rot] ee ⟩ which differ by few percent.However, we observe that for the electron flavor component, there is little variation in ⟨P [rot] ⟩ over the simulation time, and we remain close to the optically thick limit for its entirety.The differences in the solid black, solid blue, and dashed orange lines are only a few percent compared to the flavor-traced Eddington factor displayed in the dashed red curve.
In the bottom panel of Fig. 9, we again see good agreement between the black, blue, and orange curves prior to saturation.⟨P [rot] xx ⟩ encompasses a larger range than ⟨P [rot] ee ⟩.Unlike the electron flavor, here the blue and orange curves show the best agreement for the x flavor, indicating that the FLASH simulation captures the first two moments of ϱ xx in accordance with EMU over a large timespan.Furthermore, regardless of the flavor, the orange and blue curves agree closely with the black curve in the growth phase.Although we use a flavor-traced Eddington factor in FLASH, the first two moments remain accurate in that linear phase.
Fourier space analysis -Figure 10 gives the DFTs for the three NSM points at 0.1 ns before saturation, roughly corresponding to the second panel in Fig. 8.Because the growth rates are different for the second (third) points, we pick a time before saturation of 0.5 ns (0.15 ns) for FLASH, and 0.1 ns (0.05 ns) for EMU, all of which allow us to most clearly capture the fastest growing mode before it begins to non-linearly couple to other modes.For the first NSM point, we find |k| max = 6.4 cm −1 , corresponding to a wavelength of 1.0 cm and matching the distance between planar structures in the second panel of Fig. 8 (domain size L ∼ 8 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.
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  calculations for the three NSM points have larger values of Im(Ω) and |k| max compared to EMU, but still exhibit reasonable behavior characteristic of the FFI.

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 (∼ 1%) and more qualitative agreement with the last (∼ 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 DFTs are at a time prior to saturation.For NSM 1, both FLASH and EMU DFTs are 0.1 ns before saturation.For NSM 2, the DFT for FLASH occurs 0.5 ns before saturation and 0.1 ns for EMU.For NSM 3, the DFT for FLASH occurs 0.15 ns before saturation and 0.05 ns for EMU.
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).Figures 3 and 9 show that the classical MEC by itself cannot fully characterize the underlying distribution, and most certainly leads to discrepant results for the pressure tensor on the order of a few percent.As 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 implemented here requires using all of the com-ponents of E and 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 angularmoment 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 sat-uration 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;Nagakura 2023).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 generalrelativistic-magnetohydrodynamics astrophysics simulations is a major computational challenge for multimessenger-astrophysics theory.Our moment-method flavor calculations offer a contribution to this important field of study.
flavor.For r < c t the solution of the moment transport equations is that [Bottom] Errors in various components.Green and black lines correspond to difference between the calculated Eii and the theoretical predictions for Eqs.(A1) and (A2) scaled by the trace of E. Antineutrino 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 11 shows the results of this test at a time t = 3.3 × 10 −5 s.We use a mixing angle θ = 0.28818, a masssquared difference δm 2 = 6.9 × 10 −4 eV 2 , and an energy p = 1 MeV.The top panel of the figure shows the flavor content evolving with the familiar oscillatory pattern.We use the shorthand notation of E ēē ≡ E ee for the electron anti-neutrino energy density moment, and similarly for the xx 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 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, L. The dominant source of the error arises from 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 ∼ 4.1 s, the relative difference in expected versus calculated time is 1.7 × 10 −3 .At the second full period, near t ∼ 6.8 s, the drift is 2.1 × 10 −3 .Both of these drift values are within the size of the time step.
We compare this calculation using the FLASH architecture to a more direct and straightforward resolution of the equation of motion 3 , based on the study of Hannestad et al. (2006).Using the solve_ivp function from the Scipy library with high precision parameters (rtol = atol = 10 −8 ), we solve Eq. ( 12) from Hannestad et al. (2006) for the tilt angle φ.We can solve for the trace-normalized value of E ee by taking the tilt angle and solving for the Bloch polarization vector [see Eqs. ( 6) and (8) of Hannestad et al. (2006)]. 4Figure 14 gives the absolute (top panel) and relative (bottom panel) errors between FLASH and the Scipy function solve_ivp for the quantity E ee /Tr [E].The relative error is the same expression as Eq. ( A4), but we substitute the results from solve_ivp for Q (th) .The absolute error is the numerator of the rhs of Eq. (A4).We see an overall growth in the peak absolute error with increasing cycle number, consistently with the drift in the period of oscillations.Both the absolute and relative errors are the result of numerical convergence in our RKCK solver.We have verified that these errors scale linearly with the tolerance criterion of the RKCK algorithm.

B. MOMENTS IN AN ORTHONORMAL 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 (νa) , energy flux F (νa) i , the average neutrino energy ϵ (νa) , and the fluid transport velocity v i at each point in space and for 3 For the interested reader: Xiong et al. (2023) gives an analytical formulation of the curves in Fig. 13. 4 Since we consider a normal hierarchy scenario with an initial system of x (anti)neutrinos instead of electronic ones, there are some minor changes to the equations in (Hannestad et al. 2006).
Borrowing their notations:  Relative and absolute errors for Eee(t) between the FLASH calculation and a direct numerical resolution of Eq. ( 12) in Hannestad et al. (2006) .
each neutrino flavor a.The spacetime metric is given in standard 3+1 quantities: the lapse α, shift β i , and 3-metric γ 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 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)  of N ex coupled with the smaller growth rates for EMU lead to later saturation points.The second column of Table 5 gives the saturation times for all six simulations.

Figure 1 .
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 difference of blue and red vectors.[Bottom] Curves as given by Eq. (34) for νe, νe and the difference (purple) as a function of polar angle ϑ.Angular distributions for νx and νx are zero at the beginning of the simulation.Lepton number crossings occur when purple line crosses 0.
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.2for FLASH, 2-

Figure 4 .
Figure4.Magnitude of the discrete Fourier transform of Nex 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 ns before saturation, and TwoThirds 0.35 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 Table1.The light-red dashed curve has the same box-side length and half the number of grid points.

Figure 5 .
Figure 5. NSM crossing information from Foucart et al. (2016a).The snapshot is taken 5 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 } g.cm −3 .

Figure 6 .
Figure 6.Polar and cartesian representations of the initial νe and νe distributions for the NSM points.Plotting conventions are the same as Fig. 1.Note that the z axis, corresponding to ϑ = 0, is a local coordinate chosen differently at each NSM point to be coincident with the lepton number flux direction.

Figure 7 .
Figure 7. Domain averaged components for the number density moment plotted against time measured from the point of saturation.tsat differs between calculations in the same manner as Fig. 2 (see Fig.16for the same plots using directly the simulation time for the horizontal axis).Panels and curve conventions are similar to Fig.2and simulation computational parameters are given in the last two columns of Table3.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µ.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.

Figure 9 .
Figure 9. Component parallel to flux of pressure tensor plotted against time for the NSM 1 point.Line and axes conventions are the same as Fig. 3. Top panel gives electron flavor pressure tensor, while bottom panel gives heavy-lepton x flavor.The rotation of the pressure moment is different for each flavor and at each time step.

Figure 10 .
Figure10.The magnitude of the discrete Fourier transform of Nex 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 ns before saturation.For NSM 2, the DFT for FLASH occurs 0.5 ns before saturation and 0.1 ns for EMU.For NSM 3, the DFT for FLASH occurs 0.15 ns before saturation and 0.05 ns for EMU.

TFigure 11 .
Figure 11.Results from the vacuum test of FLASH plotted against radial coordinate.[Top] Diagonal density matrix components of the energy density moment (Eii) normalized by the trace of E. Barred components denote the anti-neutrino counterparts, i.e., Eēē ≡ Eee and are scaled by the trace of E.[Bottom] Errors in various components.Green and black lines correspond to difference between the calculated Eii and the theoretical predictions for Eqs.(A1) and (A2) scaled by the trace of E. Antineutrino 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 12 .
Figure12.Results from the MSW test of FLASH plotted against radial coordinate.Notation for both panels is the same as Fig.11.Error for the Eii components corresponds to the difference between the calculated Eii and the theoretical predictions using the transition probability in Eq. (A5) for neutrinos, and a similar expression for anti-neutrinos.

Figure 13 .
Figure 13.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).

Figure 15 .Figure 16 .
Figure 15.Density matrix elements versus time for the three tests in Sec.5.1.Figure content is the same as Fig. 2 except all curves begin at simulation time t = 0.

Table 1 .
Figure 3. zz component of the pressure tensor for electronflavor neutrinos plotted against time for the Fiducial test case.The pressure tensor component is normalized by the energy density moment.The solid black curve gives the pressure tensor for the EMU simulation, while the dashed red curve

Table 1 .
The light-red dashed curve has the same box-side length and half the number of grid points.

Table 4 .
Table 4 are presented in the same way as Table 2. Like the Fiducial, 90Degree, and TwoThirds tests, the FLASH Results for FLASH and EMU (2f) calculations for the three 3D NSM simulations.Column labels are the same as Table 2.

Table 5 .
Initial values of domain-averaged off-diagonal magnitude of N and saturation times for all six FLASH and EMU (2f) FFC simulations.